Interface-resolved large eddy simulation of a field debris flow
with large stones and woods
Tomoo Fukuda1*, Shoji Fukuoka1
1 Research and Development Initiative, Chuo University, 1-13-27 Kasuga, Bunkyo-ku, Tokyo 112-8551, Japan.

* Corresponding author. Tel.: +81 3 3817-1617. E-mail address: t-fukuda@tamacc.chuo-u.ac.jp.

1.

Introduction

In August 2014, an impermeable check dam was destroyed by a debris flow that occurred in Hiroshima, Japan(see Fig. 1). A large amount of huge stones and driftwood attacked a downstream residential area, causing significant
damages. In order to countermeasure against the debris flow, it is important to clarify the debris flow and its impact
force acting on the check dam. Furthermore, even when the check dam was not filled with stones, there is a possibility
that only the driftwood flows out together with the water overflowing the check dam. Therefore, it is necessary tounderstand not only motions of stones but also that of driftwood. However, it is difficult to measure mechanicalproperty of the debris flow in the field. Consequently, the impact force and the like could not be measured. In thesimulation of the debris flow, the mixture of the water and sediment is dealt in general as a continuum. In suchmodels, it is difficult to estimate the segregation of large and small stones, the impact forces of large stones onstructures, and the motion of driftwoods different from that of stones.
In our previous study [1], a numerical simulation code for the interface-resolved large eddy simulation (IRLES)
was developed. The fluid motion was simulated using the Eulerian approach with smaller computational grids thanparticles. Motions of particles of various-shapes were simulated as rigid bodies using the Lagrangian approach. For
the validation of the model, hyperconcentrated flow in a laboratory experiment was simulated [2]. The particle
motions of the simulation result agreed well with those measured in the experiment. Thus, we believe that theimplementation of IRLES of field debris flows will reveal dynamics of field debris flows and contribute tounderstand the debris flow dynamics. However, no effective parallel computation method for the IRLES has beendeveloped so far.
In this study, at first, a wood model is developed by bonding small spheres which is able to drift in the water. Thena new parallel computation technique is developed that is able to reduce computation load by automatically balancing
loads of processing elements and is effective for the IRLES of debris flows in complicated boundary geometry.
Finally, we introduced the developed wood model and the parallel computational technique to our IRLES code [2],and implemented large-scale IRLES in the upstream 200 m section of the check dam to examine the dynamics of thefield debris with large stones and driftwoods (see Fig.1).
IRLES section

2. Numerical methods
2.1. Numerical method of fluid motions
Ushijima et al. [3] developed a single fluid model for the solidliquid multiphase flows based on the IRDNS (IR direct
numerical simulation). In this study, the Smagolinsky model wasused as a sub-grid scale turbulence model for the simulation as
IRLES (IR large eddy simulation):0

Destroyed
Check Dam(1)2(2)

Figure 1. Full view of the objective valley after the debrisflow.

∆2(3)

are the volume average density and the mass average
whereandvelocity, respectively, taking into account the volume fraction of theis theparticles in a computation grid. is the kinematic viscosity,
gravitational acceleration,is the sum of the isotropic components ofpressure and SGS stress,is the strain rate tensor,
is the turbulent
viscosity, andis the Smagolinsky constant (0.1). ∆ is the computation
grid size (0.1 m). Simulations of the above equations were performed using
the finite difference method of a simplified marker and cell (SMAC) scheme
on a staggered grid arrangement. The free water surface was simulated using
the volume of fluid (VOF) method [4].

SpringFigure. 2 Wood model

2.2. Numerical method of gravel motions
The gravel-shaped particles were modeled by superimposing several small spheres without gaps. The motion of
the particles is simulated using the translational and rotational equations of motion of the rigid body:(4)(5)
where, where the index denots the components of the global coordinate system fixed on the space, the indices –denote the components of the local coordinate system fixed on a rigid body, M is the mass of a particle,
isthe position of the center of gravity of the rigid body,is the angular velocity,is the force acting on the surfaceof the particle, and Ti is the torque on the center of the gravity generated by the force. The superscripts and areis the gravitational
the components of the fluid forces and contact forces between particles, respectively,
is the inverse of the matrix consisting of the components of the moment of inertia tensor in theacceleration,
local coordinate system,

is equal to∙(andare unit basis vectors), andis the Levi–Civitasymbol.

2.3. Wood model
It is known that the bending deformation of wood makes them easier to accumulate on structures. In this study, tosimulate wood deformation, a wood model was developed by using several small spheres that were connected eachother using virtual springs (See Fig. 2).

2.4. Parallel computation method
A new parallel computation method DBA-BC (dynamic block assignment method taking into account balancing
block computation load and shortening communication times) was developed in this study. As shown in the Fig. 3,
the entire computation domain is divided into cubic domain blocks. Each MPI Processing Element (PE) computes
one or more domain blocks. At regular intervals, blocks are virtually assigned to a process in trial a plurality of times
by using random numbers. Then, from the virtual trial assignments, only a new assignment of blocks to the PEs is(see Fig. 4) by balancing computation
selected that was predicted to shorten the total computation time
among the PEs and shortening communication timethat was predicted by using counted
loads of the blocks(see Fig. 5). In this manner, this parallel computation techniquethe number of communication patterns
automatically exchanges blocks among PEs and is able to shorten computation time.

3. IRLES of the field debris flow
3.1. Simulation condition
The domain of the IRLES are shown in Fig. 6. The simulation domain was set using 1,978 cubic domain blocksof 5 m each side. The particle size distribution of the simulation was set according to the measured field gravel size′(ib,jb, kb)
jbib
Position of a block

MPI

Processing
Element

PE 1 ・・・ PE nPE 0

Dynamic block assignment to balance computationload among each processing element.

Figure 3. Schematic diagram of DBA.

: Time required for the communication
after the synchronization.
Computation time↔

kb

Boxes for storing
block data

1

Cubic computational domain block

Computation timeafter the debris flow (see Fig. 7). Particles of 5different sizes (0.4, 0.5, 0.6, 0.8, and 1.2 m) wereused in the simulation. The smallest particle sizeof 0.4 m was selected to reduce the computationalload. Using these particles, particles were packed
in to the valley to reproduce the topography before
the debris flow (see Fig. 8).
In order to determine the debris volume fromthe upstream end boundary, a two-dimensional
debris flow simulation was carried out. The flow
and the volume ratio
rate of the debris flowof solid phaseat the upstream end of theIRLES domain were obtained. However, theobtained from the twovolume ratio of solid
dimensional simulation implicitly includes thefraction of small particles like sand and silt thatwere not simulated by the IRLES. In the
preliminary IRLES in which all of volume ratio aslarge particles of 30 cm or more, the simulated
debris flow tended to accumulate in the vicinity of
the upstream end due to great resistance of thelarge particles. For this reason, the fluid phase in
the IRLES was simulated as a mixture of small
1,600kg/m .
particles and water of density
The volume ratio of stones of 0.4 m or more ,which was simulated in this study using theLagrangian approach, was determined so that thetotal volume and mass were consistent with theresults of the two-dimensional analysis as follows
1,000 kg/m ,2,650 kg/m ):(Time shortening

Ts: Time required forthe synchronization.
PE 0 PE 1 PE 2

PE 0 PE 1PE 2(a) Before the dynamic load balance. (b) After the dynamic load balance.
Rectangles indicate computational times.
Different colors indicate different blocks.

Figure 4. Schematic diagram of computation time
before and after dynamic load balance.
A number of communication patterns (Twice the number of arrows)
Equation predicting:
andare determined by the least squares method.
56

14PE2PE2PE1PE3PE 0PE 4PE5PE7PE1PE3PE 0PE 4PE5PE7PE6PE611′′′(6)(a) Blocks of the same processingelements are not adjoining.

(b) Blocks of the same processing
elsemnts are adjoining.

Fig. 9 shows the temporal changes in the flow rate Figure 5. Communication patterns depending on
the assignments of blocks to processes.
of small sediment laden flow and that of stones
supplied from the upstream end of the IRLES.

3.2. Flow configuration
Fig. 10 shows a snapshot of the front of thesimulated debris flow. By dealing small particles
and water as a mixture continuum, the simulated
debris flow was not deposited around the upstream
end and violently flew down. Furthermore, woods
were plucked out at the frontal portion of thesimulated debris flow. This made great wood

Destroyed check dam

Figure 6. Cubic blocks determining the IRLES domain.

FieldSimulation

80

Flow rate ( m3/s)

Ccumalative percentage (%)

10060
40
2000.1

1Particle size (mm)Figure 7. Particle size distribution.

450
400
350
300
250
200
150
100
500Total
Small sediment and waterStone

0102030405060Elapsed Time (s)708090100

Figure 9. Given flow rates at the upstream end.

Fixed wood set
for the roughness
Packed gravel and stones
Movable plantedwood

Reproduced Check dam
before the destruction
Figure 8. Reproduction of the bed condition before thedebris flow.

Driftwood accumulation
at the front of the debris flow.

Figure 10. Snap shot of the front of the debris flow.

accumulation at the front of debris flow. Large stones followed the wood accumulation. This phenomenon wasrecorded by the video shooting on field debris flows. Thus we believe that this IRLES was able to reproduce theflow mechanism of field debris flows.

4.

Conclusion

In this study, a deformable wood model was developed. Moreover, a parallel computation method was also
developed that was able to dynamically balance the loads of processing elements and suitable for IRLES incomplicated boundary geometries. In this way, the IRLES of the field debris flow with large stone and woods wascarried out. The simulation were able to predict wood accumulation at the front of the debris flow followed by largestones that was seen in field debris flows. In addition, this phenomenon is difficult to be explained in debris flow
continuum model. Thus, the developed debris flow IRLES model is expected to be an effective tool to examine theimpact force acting on the check dam and the barrier effect of the check dam against sediment and driftwood.
References[1] Fukuoka, S., Fukuda, T., Uchida, T., 2014. Effects of sizes and shapes of gravel particles on sediment transports
and bed variations in a numerical movable-bed channel, Adv. in Water Res., 72, 84-96,https://doi.org/10.1016/j.advwatres.2014.05.013.
[2] Fukuda, T., Fukuoka, S., 2017. Interface-resolved large eddy simulations of hyperconcentrated flows using
spheres and gravel particles, Adv. in Water Res.,https://doi.org/10.1016/j.advwatres.2017.10.037.
[3] Ushijima, S., Fukutani, A., Makino, O., 2008. Prediction method for movements with collisions of arbitraryshaped objects in 3D free-surface flows, JSCE J B, 64(2):128-138. http://dx.doi.org/10.2208/jscejb.64.128.
[4] Hirt, C.W., Nichols, B.D., 1981. Volume of fluid (VOF) method for the dynamics of free boundaries, J. Comput.
Phys., 39:201-225. http://dx.doi.org/10.1016/0021-9991(81)90145-5.