Euler-Lagrange simulation of invert
concrete abrasion
Tomoo Fukuda and Shoji Fukuoka

Abstract
This paper develops a new sub-particle scale abrasion model for the Concrete Discrete
Element Method (CDEM). The concrete abrasion model estimates the fracture rates ofconcrete particles in the CDEM method based on the mechanical energy, which themoving particles consume by collisions. The present concrete abrasion model isincorporated into the APM (Arbitrary Particle Multiphase) code (Fukuoka et al. 2004),
which is capable of simulating three-dimensional motions of arbitrarily shaped particles
and water, based on the high-resolution Euler-Lagrange strategy. The APM code with the
CDEM method and the new abrasion model (APM-CDEMA) was used to simulate
concrete invert abrasion in a real scale abrasion experiment conducted by Fukuoka et al.
(2004). The simulation results for the abrasion speed of concrete invert agreed well withthe results of the experiment.
Keywords: sub-particle scale abrasion model, DEM, high-resolution Euler-Lagrange
simulation, concrete abrasion, energy balance.

1

Introduction

Violent flows of highly concentrated water and gravel mixtures in sediment bypass
tunnels make concrete inverts abrade severely. As the abrasion progresses, sedimentflows concentrate on concave abrasion points and further increase the scale of theabrasion at these points. Therefore, the abrasion and motion of flows and particles interact
with each other. In several studies, the abrasion volume was connected with the energy
loss of moving particles during collisions with the concrete surfaces by assuming ideal
steady flows (Ishibashi 1983). However, the assumptions of the steady flows cannot
provide proper explanations for interactions between the highly concentrated sedimentflows and the abrasion.
With respect to prediction methods for water-sediment mixture flows, several studies
demonstrated the effectiveness of high-resolution Euler-Lagrange (HEL) strategies. In
the HEL strategies, particle motions are simulated using the discrete element method(DEM) based on the Lagrangian approach, and flow motions are simulated using theEulerian approach. In these simulations, fluid motions surrounding particles are simulated
accurately by using computational grids smaller than the particles, and fluid forces arecomputed directly by the surrounding flow motions without the drag coefficient. The HEL

FP14

1

Red color particles display fracture rates γf of concrete
spheres that determine the upper limits of tensile forcesamong neighboring concrete spheres.

(a)

(b)

Figure 1: Concept of the sub-particle abrasion model.
The left side model naturally simulates concrete fracture with small concrete spheres by using natural tensile
strength of concrete. However, it requires too large computation costs for the simulation. The right side modelcannot simulate directly several small particles moving out from the concrete bed, indicated in the left side
figure, but can simulate these sub-particle scale abrasions by decreasing bonding forces of concrete spheres. The
model implements concrete abrasion simulations with the CDEM method by using relatively large particles, andreducing computation load.

method can be a powerful tool for simulation of water and sediment motions in studies
on invert concrete abrasion.
With respect to the simulations of concrete behavior, Megro and Hakuno (1989)
developed a Modified DEM (MDEM) method that is capable of simulating concrete
behavior by introducing a tensile force into the DEM method. From a microscopic view
point, the phenomenon of concrete abrasion is one of concrete fractures. Thus, theconcrete DEM (CDEM) method introducing tensile forces is presumed to be a suitablemodel for predicting concrete abrasion. However, if small spheres are used for
simulations of concrete abraded powder, too many concrete spheres are necessary for thesimulation as illustrated Fig. 1(a). This makes implementation of the simulation difficult
due to the large computational load. On the other hand, if large spheres are used, collision
forces are averaged in a large space depending on sphere sizes and, estimated on averageas smaller forces than those in the contact scale. In this case concrete is difficult to fracture
without receiving extremely large collision forces. To prevent these problems, concreteabrasion simulations in the DEM require sub-particle scale abrasion models that arecapable of simulating abrasion behavior in smaller scale than particles.
Fukuoka et al. (2014) developed the APM (Arbitrary Particle Multiphase) code based on
the HEL strategy using the finite difference method and Cartesian computation grids forthe flow simulation. This method can simulate motions of arbitrarily-shaped particles, inwhich particles of various shapes are made by superimposing several small spheres. It
was confirmed in a study (Fukuda and Fukuoka 2016) that the validation of flows andparticle motions simulated by the APM was achieved by comparing them with thosemeasured in debris flow laboratory experiments conducted by Egashira et al. (2001).
In this study, the DEM concrete model with sub-particle abrasion model was developed
and added to the APM code. This simulation code is termed APM-CDEMA (APMConcrete DEM-Abrasion). If the forces and abrasion speed are connected using a subparticle scale abrasion model, an unnatural abrasion will progress in such a way that
concrete abrasion will increase because of large contact forces from a large amount ofstatic deposited particles. Our sub-particle scale abrasion model eliminates this unnatural

FP14

2

abrasion by being based on the balance between the loss of the mechanical energy ofmoving gravels and the energy of the concrete fracture. The validation of the developedmodel was examined by applying the model to the abrasion experiment conducted in astudy by Fukuoka et al. (2004). In the experiment, large quantities of real scale gravels
and water were supplied to an invert-shape concrete open channel, and concrete abrasion
shapes were measured. From the comparisons of the simulated abrasion and that in theexperiment, the applicability of the APM-CDEMA code as a prediction method of invert
concrete abrasion was discussed.

2

Methodology of the APM

The APM code is based on the high-resolution Euler-Lagrange strategy (Fukuoka et al.
2014). Water motions were computed with a one fluid model for incompressible solidliquid multiphase flow assuming solid parts as fluids of different densities and using
computational grid smaller than the particles. Gravels were made by superimposing
several small spheres without leaving any gaps. Contact forces among gravels werecomputed by the DEM with small spheres composing gravel by using springs, dashpots
and frictional sliders. Motions of gravels were simulated as rigid bodies. For more details,
refer to Fukuoka et al. (2014).

3

Concrete abrasion model

3.1 Relations between concrete abrasion and fracture in the abrasion
model
The sub-particle scale abrasion model takes into account the energy balance between theloss of the mechanical energy of moving gravels and the energy consumed to fracture
concrete in the DEM method. Ishibashi (1983) linked the abrasion volume and loss of themechanical energy of moving particles with an abrasion coefficient, as follows:Vw  Cw Εnm

[1]

where ∆Vw is the abrasion volume, Cw is the abrasion coefficient, and ∆Enm is the loss ofthe mechanical energy of moving particles. The values of the coefficients were
determined by laboratory experiments under relatively steady and simple flow conditions.
Thus, the value of the abrasion coefficient is presumed to be reliable. Therefore, if theenergy loss of the moving particles ∆Enm is predicted accurately under various conditionsof abrasions and flows, the progress of abrasion is presumed to be predicted properly.
Thus, the relation in Eq. [1] was introduced into the sub-particle scale abrasion model ofthe concrete DEM method. The concept of our sub-particle scale abrasion model isillustrated in Fig. 1(b). Eq. [1] was employed for the determination of the fracture rate of
the concrete. We will discuss this method in more detail in a later section.

FP14

3



Strength degradation
depending on fracture rates

 

Tangential force

Concrete

Gravel

Normal force(Compression)

Figure 2: Yield criterion of the DEM.

The yield criterions of normal forces and tangential forces in the DEM are shown in Fig.
2. The objective of this study is not to investigate in detail the yield criterions of concrete.
Thus, a simple, linear concrete yield criterion was used, indicated by the red solid line in
Fig. 2. The strength degradation of concrete was computed using the fracture rate andmodeled by decreasing the upper limit of the tensile force with keeping the friction angle.
The upper limit for the tensile strength of the spheres composing the concrete was variedby the following equation:)
 t ,max  (1   f ) t(,0max[2]where σ (0)t,max and σt,max are the tensile strengths at the initial time and at the correspondingtime, respectively. γf is the fracture rate. This tensile strength and the upper limit of thetensile force fb in the DEM method depicted in Fig. 2 were linked by defining the contact
area A of two concrete spheres, i and j, as follows:
fb  A t ,max, i   t ,max, j2

A  4rc2[3][4]

where σt,max,i and σt,max,j are the tensile strengths of concrete spheres i and j, respectively;
considering the strength degradations, and rc is the radius of the concrete spheres.

3.2 Fracture rate
Concrete fracture energy, which is necessary for fracturing all parts of any concrete object,is presumed to depend on the volume of the object. Thus, we defined the energy thatfractures a unit volume of concrete as CEf. Based on this relation, it is presumed thatcomplete fracture has occurred when the sum of the energy consumed in the concrete
fracture ∆Ef of volume VΩ equals to CEf VΩ. Thus, the fracture rate γf was defined asfollows.

f  Ε f
CE f V[5]

The energy consumed in the concrete fracture is assumed to be the same as the energy
consumed for abrasion of the same volume. This assumption allows connecting the

FP14

4

∆

∆

∆

∆
∆
∆
∆
∆

:Energy relative to a collision:Mechanical:Non-mechanical energy

∆
∆
∆

∆
∆
∆

:Strain energy:Energy of sound or heat
:Energy of fracture

Figure 3: Transfer of the energy relative to a collision.

abrasion volume ∆Vw and fracture energy ∆Ef by using the fracture energy per unit volume
CEf as follows.

C E f Vw  Ε f

[6]

The transfer of the energy relative to a collision is depicted in Fig. 3. The non-mechanicalenergy ∆Enm, which is the loss of the mechanical energy of the particles after collisions,is not the same as the fracture energy ∆Ef, and partially transfers into an energy ∆Esh ofsound and heat. In this energy transfer, we assumed that the ratio of the energy consumed
in fractures to the sum of the energy consumed in fractures and the energy transferring
into sound and heat was constant, and this ratio is defined αf,nm as follows.

 f ,nm Ε fΕ f  Εsh



Ε f
Εnm

[7]

The derivation of the abrasion coefficient in Eq. [1] did not consider the variation in αf,nm.
Thus, the assumption of Eq. [7] is the same as the assumption for the derivation of theabrasion coefficient. By substituting Eq. [6] into Eq. [7], the following equation is
obtained.
Vw 

 f ,nm
CE f

Εnm

[8]

By comparing Eq. [1] and Eq. [8], the following relation between the fracture energy CEfof a unit volume and the abrasion coefficient Cw is obtained.
Cw  f ,nm[9]CE f

By substituting Eq. [7] and Eq. [9] into the definition of the fracture rate in Eq. [5], thefollowing equation is obtained.

f 

 f , nm  Εnm
C E f V

C w  Εnm
V

[10]

Thus, when the loss of the mechanical energy of particles ∆Enm is known, the fracture rateγf can be computed with the abrasion coefficient Cw.
The energy, transfer from the mechanical energy of contacting two particles i and j into
non-mechanical energy in collisions, can be obtained by the following equation:

FP14

5

Hopper for
sediment suply

Invert shapeconcrete open
channle(a) Full view of the experimental setup

Flow(c) Abraded concrete invert after
3,500m3 sediment was supplied(b) Sediment supplying section.

Figure 4: Photos of the invert concrete abrasion experiment.
z

48m
45m
−2.75m−1.75m−0.5m 0.5m
Water supply
Sediment supply
sectionsectionComputatin grid : 0.01m

y
x

x
45m

Figure 5: Numerical setup.
1.06m

Figure 6: Arrangement of concrete spheres in a cross section.

 Ε nm ,ij  f c ,i  ( v j  v i )  t[11]where fc, i is the contact force vector from the particle j to the particle i, vi and vj are thevelocity vectors of contacting points of the particles i and j. When the material strainenergy of the particles is converted into mechanical energy, ∆Enm,ij becomes a negativevalue. However, mechanical energy is not produced during one total collision. Therefore,
the sum of ∆Enm,ij becomes a positive value and is able to be used for predicting the lossof the mechanical energy of particles. In the simulation, this ∆Enm,ij was divided withrespect to particles i and j as ∆Enm = 0.5∆Enm,ij, respectively.

4

Simulation of invert concrete abrasion experiment

4.1 Numerical setup
Fukuoka et al. (2004) conducted an experiment on invert concrete abrasion. In theexperiment, a 45 m long invert-shape concrete open channel with a bed slope of 1:20 wasused. A large quantity of gravels and water were supplied from the upstream, and changing
shapes of the invert bottom due to the abrasion were measured. Figure 4 shows a full viewof the experimental set up (left), a view of the sediment supplying section (middle) and theabraded invert bottom (right). The developed APM-CDEMA code was used in thesimulation of this invert concrete abrasion experiment. From a comparison between the

FP14

6

Accumalative percentage
(%)

100

Exp.
Cal.

80
60
40
20
0
0.01

0.1

1
10
Particle size (mm)

100

1000

Figure 7: Particle size distributions of the experiment and the simulation.

Red color of the bottom indicats
that fracture rate is great.

Figure 8: Snap shot of gravels, water surface and invert abrasion.

progress of concrete abrasion in the simulation and that in the experiment, the validation of
the result of the APM-CDEMA was examined. The numerical set up was shown in Fig. 5.
The cross sectional view of the numerical concrete invert is shown in Fig. 6. The
experimental concrete invert was represented by groups of spheres of 0.02 m in thesimulation. The particle size distributions of the simulation and that of the experiment wereshown in Fig. 7. Particles with 5 sizes (0.04 m, 0.05 m, 0.075 m, 0.1 m, and 0.15 m) wereused in the simulation. The rate of water supply was 0.5 m3/s and that of the sediment was
0.0083 m3/s, same as in the experiment. Although, the experiment was performed for 80
hours, we could only simulate the first 33.3 hours in the experiment because of large
computational load. Moreover, we shortened 33.3 hours of the experimental time to 25 s ofthe simulation time. To compare the simulated abrasion with that of the experiment, thenon-mechanical energy ∆Enm in Eq. [11] was multiplied by the acceleration parameter β (=4,800(≃33.3×60×60/25)). The abrasion coefficient Cw = 1.1×10−7 m3/J and the initial
tensile strength σ (0) t,max = 2.7×106 (Pa), which was 1/12 the compressive strength of theexperimental invert-concrete, were used in the simulation.

4.2 Comparison between the abrasion in the simulation and those in theexperiment
A Snap shot of gravel, water surface and invert concrete abrasion in the simulation isdepicted in Fig. 8. Fracture rate of the concrete is shown as red color. Large fracture rates
prevailed in a wide range of the channel bottom, which implied that sub-particle scaleabrasion model has great effect on the abrasion simulation. Because of the short

FP14

7

y (m)

Abrasion
depth
(m)
x (m)

y (m)

Flow

x (m)

* The abrasion depth in the simulation was measured by determining radii of the abraded concrete spheresas (1− f )rc.

Figure 9: Contour of the abrasion depth

Abrasion Speed
10-5 (1/m2)(Top: Experiment (At the time of the sediment volume of 245 m3 supplied), Bottom: Simulation (At thetime of the converted sediment of 91 m3 having passed the section x = 10 m))

4
3
2
1
0
Exp.

Cal.

Figure 10: Abrasion speed at the center of the channel.

simulation time of 25 s, at the end of the simulation, many particles still remained at the
supplying section. Therefore, we compared concrete abrasion of the simulation and thatof the experiment in the relatively upstream section from x = 10 m to 20 m. The abrasiondepth contour of the experiment (at the time where the supplied sediment volume was
245 m3, which is the first measurement after the sediment supply) and the simulation (atthe time of the converted sediment of 91 m3 passing the x = 10 m section) are shown in
Fig. 9. The sediment volume of the simulation converted to the experimental volume of
91 m3 was determined by multiplying the simulated sediment volume passing the section
x = 10 m by the acceleration parameter β. Because both channel beds of the simulationand of the experiment do not yet display wavy abrasion by this amount of sedimentvolume, motions of the water and the sediment were not affected greatly by the abradedchannel bed. A comparison of the speed of abrasion depth per unit passing sedimentvolume at the center of the channel in the section from x = 10 m to x = 20 m is shown in
Fig. 10. The abrasion speed in the simulation matched well with that of the experiment.
Therefore, simulated moving particle, concentrating on the center of the channel, and therelation between the progress of the abrasion and the particle motions presumed to be agood match with those in the experiment.

5

Conclusion

In this paper, APM-CDMA code was developed, in which motions of water and particles
were simulated in the HEL strategies, and concrete was simulated in the CDEM methodwith the new sub-particle scale abrasion model. The developed APM-CDEMA code wasapplied to the experiment of the invert concrete abrasion (Fukuoka et al. 2004). Because
the simulation time was still short, the model validation could be examined with only theearly stage of the abrasion in the experiment, in which bed waves of the invert due to the

FP14

8

abrasion were not confirmed yet. The examination brought that the abrasion speed
simulated with the APM-CDEMA code agreed well with that of the experiment in the earlystage of the invert abrasion. We will perform this simulation for longer time hereafter, andfurther validations should be checked. The APM-CDEMA code is capable of simulatingcomplicated interactions between motions of sediment flows and the abrasion. This HEL
strategy is expected to be a powerful tool for predicting the invert concrete abrasion.

References
Ishibashi, T. (1983). A hydraulic study on protection for erosion of sediment flush equipments of dams,
Proc. JSCE, 334(6), 103–112 (in Japanese).
Egashira S., Itoh T., Takeuchi H. (2001). Transition mechanism of debris flows over rigid bed to overerodible bed, Phys.Chem. Earth (B), 26(2), 169-174.
Megro, K., Hakuno, M. (1989). Fracture analyses of concrete structures by the modified distinct element
method. Proc. JSCE, 410, 113-124.
Fukuoka, S., Shinohara, Y., Masaki, J. (2004). Shigemura, K., Todo, M., Okada, S., Saito, K. : Experimental
study on the development mechanism of abrasion in a concrete channel invert due to flushingsediment-laden flow., Proc. Hydraul. Eng., 48, 1135-1140 (In Japanese).
Fukuoka, S., Fukuda, T., Uchida, T. (2014). Effects of sizes and shapes of gravel particles on sedimenttransports and bed variations in a numerical movable-bed channel, Adv. Water Resour., 72, 84-96.
Fukuda, T., Fukuoka S. (2016). Euler-Lagrange simulations of debris flow experiments, Proc. Two-pHase
modElling for Sediment dynamIcS in geophysical flows THESIS-2016, Tokyo – Japan, 17-20.

Authors
Tomoo Fukuda (corresponding Author)
Shoji Fukuoka
Research and Development Initiative, Chuo University, 1-13-27 Kasuga, Bunkyo-ku,
Tokyo 112-8551, Japan
Email: t-fukuda@tamacc.chuo-u.ac.jp

FP149