Three-dimensional simulation of non-uniform sediment transport based on multi-phase Eulerian approach: Application to debris flow

To simulate 3D flow of a non-uniform and highly concentrated sediment, a numerical model using a multi-phase Eulerian method for air, water, and particles of various class size is developed. This model accounts for turbulence in pore water, particle-particle collisions, and enduring friction. To test its performance, simulations were performed for large scale debris-flow experiments. The numerical model reproduces successfully the flow property and sediment size segregation of the debris flow.


Introduction
In the last few decades, three-dimensional (3D) numerical simulations of solid-liquid multi-phase flow has been actively discussed. However, difficulties remain in simulating highly concentrated sediment-laden flow such as debris flow, because such flows cannot be represented by means of a bed-material load (i.e., bed load or suspended load). Recently, the discrete element method (DEM) coupled with a fluid dynamics model was successfully applied to simulate such flows. For example, Fukuda et al. [1] simulated with some success a 3D debris flow by coupling DEM, large eddy simulation (LES), and volume-of-fluid (VOF) method, albeit at a heavy computational cost. A more practical numerical method is required to model an actual flow field. The Eulerian-Eulerian approach is also used in 3D simulations of sediment-laden flow such as sheet flow [1] or local scour around a hydraulic structure [3]. These models use momentum equations for both water phase and particle phase, including the closure of water-particle and particle-particle interactions. As the Eulerian-Eulerian approach needs to define no bed surface, it has the advantage in that the water-sediment interface is ambiguous (e.g., sheet flow). Furthermore, the Eulerian-Eulerian approach enables the flow velocities of water and particles to be computed separately. Unlike in a typical one-fluid model, each dynamics of water and particles is resolved, which offers advantage when water and sediment behave differently.
However, previous studies on sediment transport employing the Eulerian-Eulerian approach assume the sediment is of uniform size. Hence, sediment size segregation is not simulated in these studies. The sediment size segregation in debris flow strongly affects its flow properties and collisional forces on structures [1]. From a practical scientific perspective, the numerical model needs to be adapted to simulate non-uniform and highly concentrated sediment transport. Our study has therefore focused on developing such a 3D multi-phase Eulerian model for non-uniform and highly concentrated sediment-laden flow. The model was validated against several analyses of experiments involving large-scale debris flows.

Numerical model
The basic equations consist of the continuity equation and momentum equation for the water phase and the particle phase, the air phase being considered only in the continuity law: (5) where subscript w signifies the water phase, k (= 1, 2, …, N) indexes the N particle phases, α denotes the volume fraction (α air the volumetric fraction of the air phase), c s the particle concentration given as a sum of volume fractions of particle phases, ∇ the 3D gradient operator; U velocity vector; p the pressure; g the gravitational acceleration; S the strain rate tensor; the coefficient representing momentum exchange between water phase and particle phase. ρ density, the stress tensor of water phase; and are the collisional and frictional stress between particles, respectively, both of which are distributed among the momentum equations for the particle phase as determined by the volume fraction of the particle-size class α s,k .
When the VOF is applied to the continuity equations [Eqs. (1) and (2)], α becomes known and the velocity U and pressure p remain as unknowns in each momentum equation. As the pressure p commonly acts on the water and particle phases, the sums in Eqs.
where is the viscosity of water, the eddy viscosity for the water phase, the local average diameter in the computational cell, and D the strain velocity tensor. The stress tensors concerning the particle phase in Eq. (4) include the collisional stress and the frictional stress . The form of the collisional stress comes from [5] and is based on the Boltzmann equation. The present study considers the collisional pressure S produced in the collision between two particles, which have a local average diameter and a particle concentration c s . Using the collisional pressure and the strain velocity tensor, the collisional shear stress for every particle-size class is computed from where tr denotes the trace of tensor, dev the deviation tensor, I the unit matrix, the collisional pressure, the quasi-viscosity generated by the particle collision, the volumetric viscosity, e the coefficient of restitution, T the quasi-particle temperature, which is related to the fluctuating energy of the particles with size and is usually computed based on their transport equation. As the focus is on the highly concentrated flow, we assume a local equilibrium for the source and sink terms, , (11) , (12) where the dissipation rate of T, the particle velocity weighed and averaged by the volumetric fraction of particle-size class, the strain velocity tensor using , the collisional stress computed using (see Eq. (8) for the collisional stress), and g 0 the radial distribution function, which is function of c s [8]. Rearranging Eqs. (9)-(12) yields an algebraic equation in T; for exact details see [6].
The frictional stress represents the Coulomb friction of the massive particles and emerges in the highly concentrated flow. The frictional pressure (i.e., normal stress) is calculated using the particle concentration c s and local average particle diameter . The frictional shear stress is computed using the strain velocity for each particle-size class.
Considering the dilatancy of soil, is given as [7][9], , , (14) , (15) where is the frictional pressure, the quasi-viscosity due to friction, the frictional coefficient between particles, and I 0 a constant having value 0.28. When the strong shear stress acts on the high concentrated massive particles, these massive particles dilate and the inertial frictional coefficient increases. This dilatancy is represented in Eq. (14). When the dimensionless inertial number I s increases, also increases from a minimum value to a maximum value . The present study adopts and based   [10]. During dilations, the particle concentration c s decreases with I s , , (16) where c s,min the particle concentration at which the dilation begins to occur, c s,max the maximum particle concentration, and b a constant. Substituting I s of Eq. (15) into Eq. (16) yields . (17) Considering previous works [7], the present study sets c s,min =0.5, c s,max =0.65, and b=2. Note that, by using , Eq. (13) can be converted to the stress tensor of the Coulomb friction, .
The final terms on the right-hand side of Eqs. (3) and (4) physically represent the drag between water and particle. The coefficient is computed using the well-known formula given by Schiller and Naumann [11], where C D the coefficient of drag for the particle; for details, see [11].

Results and discussion
The numerical model was validated in a simulation of a large-scale debris-flow experiments, from which we had extracted experimental data of the flow, deposition, and sediment size segregation [11] [13]. The experiment was conducted using a flume 95 m-long, 2 m wide, and a bed slope of 31°. For the initial conditions, a non-uniform sediment with d 50 =5.0 mm was released upstream from the flume. The debris flow is deposited in a run-out (bed slope of 4°) extending from the lower end of the flume. Two experimental conditions were conducted: one flume with a rough surface (16-mm-height and 50-mm-space) and the other without rough surface.
In the simulation, a hexahedral mesh were set parallel to the lower bed. The mesh size was 20 mm in the normal direction and 100 mm in the tangential direction to the lower bed. The particle diameter distribution was divided into five classes. Initially, a uniform concentrated and saturated sediment of 10 m 3 was set at upstream end of the flume as well as the experiment. The experiment and simulation were started as dam-break flow. Non-slip conditions were imposed on the wall boundary of the water flow velocity. Regarding the wall boundary of the particle flow velocity, two different conditions were tested: a non-slip condition, and a partial-slip condition as stipulated by Johnson and Jackson [14]. They assumed a slip velocity to balance the tangential momentum in particle-wall collisions, specifically , (20) where x n is the axis normal to the wall, and the coefficient of specularity representing the rate of momentum change a particle-wall collision and takes a value within . The slip velocity decreases with increasing . Takahashi and Tsujimoto [15] used this partialslip condition coupled with a one-fluid Eulerian model, to examine the bed material and .
With this insight, the present study used =0.20 for a concrete bed, without any roughsurface element, and was used in several runs of the experiment. This boundary condition was an intent to take into account the effect of roughness which is far smaller than the scale of a computational cell.
In the experiment, the front of the debris flow was tracked using image analysis from video clips. Fig. 1 shows temporal variation of the front location of the debris flow. The experimental results are from two runs with and without rough-surface element. Each line represents an average computation from data obtained from conducting ten experimental runs. The simulation results of the front location represent the location where the sum of the volume fractions for water and sediment reached 0.5. Obviously, the simulation under the non-slip condition shows similar behaviour as the experiment with roughness. Similarly, the simulation under the partial-slip condition shows similar behaviour as the experiment without roughness. For both experiment and simulation, the difference between runs with/without roughness diverges after x = 30 m. With roughness, the flow decelerates after 30 m, whereas with no roughness the velocity remains nearly constant.
Iverson et al. [13] reports on experiments run with rough-surface elements: After 3-4 s, coarse grains coalesce at the flow front to form a gravel snout as a result of migration of surface grains toward the flow front. Conversely, the debris flow trailing the snouts appeared to be distinctly finer and wetter. When debris-flow snouts were observable, however, they exhibited wave-breaking forms in which clasts reaching the crests of snouts tumbled down their forward faces, slowed when they contacted the bed, and were overridden. This dynamic behaviour of the debris flow was well captured by the numerical simulation. Figure 2 shows the local average diameter of the flow in the longitudinal plane. For clarity in visualization, the longitudinal plane is rotated horizontally and amplified in the direction normal to the bed. The rigid-line enclosing the debris marks αair = 0.5. In the first ~4 s, coarse grains coalesce at the surface. As observed in the experiment [13], gravel snouts form at the front. During the interval 4-6 s, the front shape becomes precipitous. This is because the gravel snouts strongly resist the rough-surface element on the bed and the subsequent flow pushes the gravel snouts forward. This process causes the deceleration after ( Fig. 1). During the 6-8 s, the grave snouts mainly float above the water (i.e., outside of the rigidline showing α air = 0.5).This result corresponds to the experimental fact that the gravel snouts barely retain their water and the grains saltate in air. In addition, subsequent flow clearly moves faster than the front flow during the 6-8-s period, which also corresponds to the experimental fact that the gravel snouts are pushed forward by the subsequent flow. Although the experiment and simulation somewhat deviate after 8 s, the experiments of this temporal stage have large scatter and shows a 10-m difference between the maximum and minimum of the repeated experiments. Considering this experimental scatter, the numerical model provided a reasonably simulation of the velocity of the front flow. However, in the first 2 s, the simulation results for both non-slip and partial-slip conditions produced faster flow velocities than those in experiments. The reason probably is that suction suppresses the initial motion of the unsaturated soil. Figure 3 shows a comparison of the flow depth measured at . The experiment with the rough surface shows increases during the period 4-7 s and the simulation of the nonslip condition shows a peak at 6-7 s (see Fig. 3a). This is because of dam up of the subsequent m 32  x flow, which is caused by the deceleration of the front flow. The experiment in the absence of roughness shows a peak at 4-5 s and decreases afterwards, and is mirrored in the simulation result under partial-slip condition (see Fig. 3b). Thus, in both experiment and simulation, a dam up of the subsequent flow does not occur when the roughness is weak. Unfortunately, the simulation overestimates the velocity of the subsequent flow. This might be because the mixture length defining the turbulent length scale is assumed to be the porous scale between particles. As the trailing flow consists of very fine particles, the eddy viscosity of the pore water is underestimated as a result. Figure 4 shows the deposition after the debris flow went out of the flume. The coloured particles in the experiment are coarse pebbles (2 cm). Colour on the simulation represents the local particle diameter on the surface of the debris flow. Although the simulation shows a larger deposition area than the experiment from an overestimate of the velocity of the subsequent flow, the simulation yields a horseshoe-shaped deposition of coarse particles similar to one observed in the experiment. This size segregation comes from a spreading out of the grave snout and the deposition of finer sediment on the gravel. All of the numerical results confirm the validity of the present 3D numerical model for highly concentrated sediment flow with non-uniform particles.

Conclusion
In this study, we developed a multi-phase Eulerian model for non-uniform and highly concentrated sediment-laden flow. The model was validated against two experiments of large-scale debris flow. The numerical model well captured the velocity, flow height, and sediment size segregation of the unsteady debris flow with and without artificial roughness. Sediment size segregation occurs not only over the longitudinal plane but also transversely, especially during deposition. The 3D size segregation was also captured in the present numerical model. Although this work focused on debris flows, the present numerical model also has the advantage in accounting for the significant difference in velocity between water and particle. One example for its application is the phenomenon of sediment abrasion on hydraulics structures of hydroelectric facilities. The present model thus may have wide potential in contributing to engineering assessments related to rivers.