Proceedings of Coastal Engineering, JSCE, Vol.3, 2012, November

NUMERICAL ANALYSIS ON DEFORMATION OF
WAVE DISSIPATING BLOCKS BY
GPU-ACCELERATED DEM COMPUTATION
Naoki Yoshifuji1 , Hiroyuki Ikari2 and Hitoshi Gotoh3
1

Dept. of Civil and Earth Resources Eng., Kyoto Univ., yoshifuji.naoki.78n@st.kyoto-u.ac.jp
Dept. of Civil and Earth Resources Eng., Kyoto Univ., ikari@particle.kuciv.kyoto-u.ac.jp

23

Dept. of Civil and Earth Resources Eng., Kyoto Univ., gotoh@particle.kuciv.kyoto-u.ac.jp

In a breakwater covered with wave dissipating blocks, a leg of a wave dissipating block and wall ofcaisson can be broken due to high waves. To simulate this kind of phenomenon, a numerical model isrequired to treat deformation and failure of a block. In order to simulate them in detail, a block must
be composed of sufficient number of elements, and a time increment should be set smaller to satisfy
CFL condition which is limited by allowable tensile displacement. These requirements lead to a highcomputational load. Therefore, in this study, GPGPU with OpenCL is introduced to accelerate DEM
computation. Deformation and failure process of wave dissipating blocks under high waves is simulated.
Key Words : Wave dissipating block, Distinct Element Method, OpenCL, GPGPU

1. INTRODUCTION
A breakwater covered with wave dissipating
blocks is one of the most popular breakwaters. In
these days, some analyses on a wave dissipating
block have been executed by DEM, however, blocks
in previous studies were not always suitable for a
detailed investigations because they were composedof a few elements or treated as a rigid body.
Passively moving solid model (Koshizuka et al.,1998) can be applied to connect DEM particles totreat as a rigid body. Gotoh et al. (2005) simulatedcompaction process of wave dissipating blocks due tohigh waves by this model. However, in fact, because
failure of a leg of block or a caisson wall can occur, it
is desirable to treat deformation and failure of a block.
A block should be composed of sufficient numberof particles for accurate simulation. A time increment
should be set smaller based on CFL condition,
because a block made of concrete has very small
allowable tensile strain. This fact lead to a highcomputational load. Heterogeneous computation, in
which not only CPU but also a processor other than
CPU (e.g. GPU) are used, is expected to accelerate
DEM computation and to enable to analyze wave
dissipating blocks which can be deformed and fail.
In this study, GPGPU with OpenCL is introduced
to accelerate DEM, and deformation and failure

Fig. 1 Kelvin-Voigt model

process of wave dissipating blocks in a breakwater issimulated.

2. COMPUTATIONAL METHOD(1) Distinct Element Method
Distinct Element Method (DEM) is a numericalmodel proposed by Cundall and Strack. (1979) basedon an overlap between two elements. A force isevaluated by Kelvin-Voigt model, which is composedof elastic spring and viscous dash-pot as shown
Fig.1. A spring constant and viscous modulus are
based on Hertzian elastic collision theory and thecritical condition of 1-dimensional damping theory as
follows:
di djE2(1)δnkn =
3 1 − ν22 (di + dj )- 16 -
 












Fig. 2 Force and displacement

1kn
2 (1 + ν)cn = 2 M kn
1cncs = 2 (1 + ν)
ks =

(2)
(3)
(4)
Fig. 3 OpenCL Device

where kn is the spring constant for normal force, ks is
the spring constant for shear force, cn is the viscous
modulus for normal force, cs is the viscous modulus
for shear force, E is the Young’s elastic modulus, νis Poisson’s ratio and M is the mass of a particle.
Young’s modulus and Poisson’s ratio are not equal
to real material constants but tuning parameters ofDEM.
(2) Deformation of concrete
The model proposed by Katsuki and Haraki(2009)is applied to simulate deformation and failure ofconcrete in this study. Fig.2 shows relation between
force and displacement in this model.
Normal spring force is calculated by Popovics’s
formula as follows:N qn N
N −1+
qn max fc (qn ≥ 0)
fn max =ft (qn < 0)
qc (qn ≥ 0)
qn max =qt (qn < 0)fn = fn max

qn
qn max

(5)

(6)
(7)

where fn and qn are normal force and displacementof spring respectively, fn max and qn max arenormal strength and maximum displacement ofspring respectively, fc and ft are compressive andtensile force corresponding to strength of concrete
respectively, qc and qt are compressive and tensile
displacement corresponding to maximum strain ofconcrete respectively and N is a parameter for
non-linearity.
Shear spring force is calculated with plastic slide

by Mohr-Coulomb failure criterion as follows:
fs = min (ks qs , C + fn tan φ)(8)

C = αC0(9)

tan φ = α (tan φ)0⎧1qs ≤ qs y⎪⎪⎨(1 − αy ) qs p
α= 1− q −q
qs y < q s < q s fsfsy⎪⎪⎩qs ≥ qs fαy(10)(11)

qs p = qs − qs y(12)qs f = rγs f(13)

where fs and qs are shear force and displacement ofspring respectively, C is cohesion, tan φ is coefficientof friction (φ is friction angle), α is reducing
coefficient, C0 is initial cohesion, (tan φ)0 is initialcoefficient of friction, qs y is shear displacement
when shear force reaches the Mohr-Coulomb failure
criterion, qs f is shear displacement corresponding tomaximum shear strain of concrete, qs p is shear plasticdisplacement, r is radius of particle element, αf is
the minimum value for reducing coefficient and γs yis deformation angle after softening. This model
describes that a critical value of concrete is reducedby plastic strain.

3. GPGPU WITH OPENCL(1) OpenCL
OpenCL (Open Compute Language), which ismaintained by Khronos Group, is an open andstandard framework for parallel computing using any
processors including CPU, GPU, DSP and Cell. And
it does not depend on platform. OpenCL ver.1.1,
which is the latest version including C++ bindings,
is used in this study.

- 17 -

Table 1 Computing environment
Microsoft Windows XP Professional SP3Intel Xenon E5420
3.00[GB]NVIDIA GeForce GTX 295
ver. 6.14.0012.9610

         

Each processor (e.g. CPU or GPU) is treated as
one of the OpenCL Devices. The OpenCL Platform
is composed of OpenCL Devices and a hostprocessor operating them. Fig.3 shows a compositionof OpenCL Device.
Arithmetic processing is
executed parallelly by Processing Element (PE)in Computation Unit (CU). There are 4 kinds ofmemories in OpenCL. Each PE has Private Memory,which can be accessed only by the owner PE. Local
Memory is shared with all the PEs belonging to same
CU. Global Memory is shared with any PEs and CUs
belonging to same Device. Constant parameters are
stored in Constant Memory.
Programs for OpenCL are categorized into twogroups; Host Program and Kernel Program. HostProgram is executed on a host hardware (usually asingle core of CPU) and Kernel Program on OpenCL
Devices (usually CPU cores other than host or ahardware other than CPU). Firstly, Host Program
creates an executing environment called Context
based on information provided from Platform.
Context is composed of following objects:• Command-Que, which operates devices• Memory object, which handles data on Device• Kernel object, which means one process onKernel Program• Program object, which is composed of more thanone Kernel object
Kernel Program is written in OpenCL C language,
which is based on C99. Host Program operates Kernel
Program with OpenCL APIs.
(2) Acceleration rate of computing
Computing environment in this study is shownin Table 1. Free-fall to floor by the same 80t
block as one used in Chapter 5 is simulated tomeasure computing time. Acceleration rate (ratio ofcomputing time of CPU to GPU) in 4 cases (numberof blocks N = 1, 8, 34, 64) are shown in Fig.4.
Blocks arranged straight are fallen down to floor at
the same time. Single block is composed by around
700 number of particles.
Each calculation process shows different
acceleration rate. The more number of particles
becomes, the more efficiently acceleration rate isfound. ”ForceParticles” (which means calculation  OSCPURAMGPUGPU driver   

Fig. 4 Acceleration rate






Type
0.5t
5t
80t

d967
2074
5375

S
1075
2306
5975

h900
1930
5000

2r1
270
580
1500

2r2
540
1160
3000

Fig. 5 Overview of wave dissipating blocksof inter-particle forces) and ”MoveParticles” (whichmeans update of positions and velocities of particles)
increase remarkably in comparison to ”UpdateGrid”(which means to store particle index to computational
cell on the particle location) and ”DetermineDt”(which means determination of time increment Δt tosatisfy CFL condition). This is because the former
processes depend on the number of particles morestrongly than the latter processes. The maximum
accelerations is about 35 times on ”UpdateGrid”
in case of about 1 million particles. However,accelerations of total computing time is to about
25 times in 1 million order particles case because”ForceParticle” process occupies most of total
computation time. Therefore, it is most important toaccelerate this process.
Double-precision floating point calculation,
computation speed of which is expected slower thansingle-precision floating point, is applied in thisstudy to simulate deformation of concrete, allowable- 18 -    

 

Fig. 8 Horizontal collision and failure (black line)
Case I: v = 1.2 [m/s], Case II: v = 3.0 [m/s] 





Fig. 6 Collisional velocity v and collisional force F

Fig. 9 A breakwater after packing computation
Fig. 7 Concrete plate[mm]displacement of which is very small. This is thereason why acceleration rate is not so higher incomparison to result by Harada et al. (2007) although
3 generation newer GPU than their one is used in thisstudy.

4. COLLISION OF SINGLE BLOCK(1) Collisional force to rigid wall
To determine Young’s modulus in Hertzian elastic
collision theory, a wave dissipating block (0.5ttype shown in Fig.5) colliding with a rigid wall
is simulated. Collisional velocity and force areshown in Fig.6. Calculated results in case of E =90 [GPa] show good agreements with experimental
result (approximate lines by least-squares method)
by Arikawa et al. (2005). Therefore, this value isadapted in this study.
(2) Failure of concrete plate
To examine applicability of the method proposedby Katsuki and Haraki (2009), a 5t type block

colliding with a concrete plate is simulated. The
concrete plate is shows in Fig.7. This is almost
same as experimental one by Arikawa et al. (2005),
however, reinforcing steel bars are located on only thecenter of plate and treated as a rigid body to acceleratecomputation speed.
Fig.8 shows results when a block collides witha concrete plate horizontally. In lower speed case,vertical cracks are found on the back face of theconcrete plate.
This means the occurrence ofbending failure. On the other hand, in higher speedcase, there are many failures intensively around thecollision point. In this case, punching shear failure
occurs. These results are in good agreements withexperimental ones by Arikawa et al. (2005).

5. ANALYSIS ON BREAKWATER UNDERHIGH WAVES
Fig.9 shows an overview of the breakwater
covered with 80t type wave dissipating blocks after
packing computation. Wave force is given to blocks
as an external force changing periodically (1.0 [s]loading and 4.0 [s] unloading, wave height H is- 19 -Fig. 11 Rocking motion and failure of caisson wall
(wall is shown as the dotted plane, block-colored onemeans failure)Fig. 10 Failure of a tip of block(black-colored particle means occurrence of failure)10.0 [m] and water depth in front of breakwater h
is 15.0 [m]). Wave force profile in vertical direction
is based on the Goda’s formula. These conditions
are the same as Gotoh et al. (2005). Width andheight of the breakwater are also almost the same,
however, the location of each block is different
because of dependency on initial condition beforepacking. Therefore, 5 cases, in which initial positionof each block is changed randomly in a scale as half
length of leg of block, are executed.
Fig.10 shows a block on the top layer drops down
and collides with another block. Fig.11 shows a blockrocks and collides with caisson wall.

6. CONCLUDING REMARKSIn this study, DEM was accelerated by GPGPU fordeformation and failure of concrete. The maximum
acceleration rate reached 25 times. Results ofcollision of a block showed the same failure type as
experimental one. In a simulation on a breakwatercovered with wave dissipating blocks under highwaves, failures which caused by dropping down androcking of block were found.
In the present model, it is possible to executestable calculation continuously even if a block splits
into more than two bodies, however, any splits did
not occur in this study. This would be because
complicated flow in void of blocks was not treated.

In a future, effectiveness and advantage of thepresent model would be shown by coupling with afluid dynamics model.
REFERENCES
Cundall, P.A. and Strack, O.D.L (1979): A discrete
numerical
modelforgranularassemblies,Geothechnique, Vol.29，No.1，pp.47-65.
Koshizuka, S., Nobe, A. and Oka, Y. (1998): Numericalanalysis of breaking waves using the moving
particle semi-implicit method International Journal for
Numerical Methods in Fluids, Vol.26.
Gotoh, H., Harada, E., Takayama, T., Mizutani, M.,
Fudou, M. and Iwamoto, T. (2005): Compaction
mechanism of wave dissipating blocks due to highwaves, Annual Journal of Coastal Engineering JSCE,
Vol.52, pp.781-785. (in Japanese)
Katsuki, S. and Haraki, D., (2009): Simulation of fragment
scattering of concrete specimen subjected to high speed
loading by Distinct Element Method, AIP Conference
Proceedings, Vol.1145, Issue.1, pp.751-754.
Harada, T., Tanaka, M., Koshizuka, S. and Kawaguchi,
Y. (2007): Acceleration of rigid body simulation
using graphics hardware，Proceedings of the 2007
symposium in Interactive 3D Graphics and Games.
Arikawa, T，Ikebe, M.，Osaki, N., Kuroda, T., Orita,
T. and Simosako, K. (2005): Study on punching
shear failures of caisson wall caused by the impulsive
force of dissipating blocks Report of the Port andAirport Research Institute，Vol.44，No.1，pp.43-83.
(in Japanese)- 20 -(Received June 15, 2012)