**Numerical Simulation of Non-Homogeneous Viscous Debris-Flows Based on the Smoothed Particle Hydrodynamics (SPH) Method**

#### **Shu Wang 1, Anping Shu 1,\*, Matteo Rubinato 2, Mengyao Wang <sup>1</sup> and Jiping Qin <sup>1</sup>**


Received: 11 September 2019; Accepted: 1 November 2019; Published: 5 November 2019

**Abstract:** Non-homogeneous viscous debris flows are characterized by high density, impact force and destructiveness, and the complexity of the materials they are made of. This has always made these flows challenging to simulate numerically, and to reproduce experimentally debris flow processes. In this study, the formation-movement process of non-homogeneous debris flow under three different soil configurations was simulated numerically by modifying the formulation of collision, friction, and yield stresses for the existing Smoothed Particle Hydrodynamics (SPH) method. The results obtained by applying this modification to the SPH model clearly demonstrated that the configuration where fine and coarse particles are fully mixed, with no specific layering, produces more fluctuations and instability of the debris flow. The kinetic and potential energies of the fluctuating particles calculated for each scenario have been shown to be affected by the water content by focusing on small local areas. Therefore, this study provides a better understanding and new insights regarding intermittent debris flows, and explains the impact of the water content on their formation and movement processes.

**Keywords:** non-homogeneous debris flow; viscous coefficients; intermittent debris flows; energy conversion

#### **1. Introduction**

As a frequent natural disaster, debris flow poses a serious threat to the lives and properties of people living in mountain areas [1–4]. This natural phenomenon needs to be taken into consideration when planning new developments in mountain catchments because its sudden, devastating, and extensive impacts could have strong consequences on the local economy. Debris flows can be divided into multiple categories such as viscous flow, dilatant flow, dilute flow, water-rock flow, etc. [5] based on the variety of materials and their combinations. The first category mentioned, the viscous debris flow, is particularly widely distributed in the mountainous areas of the southwest of China. This particular debris flow is characterized by high viscosity, wide particle-size distributions and uneven velocity distributions [6]. The difficulty in describing its motion process is due to the fact that the viscous debris flow is classified as a heterogeneous, non-constant, and non-Newtonian flow [7]. In previous studies, Johnson [8] attempted to describe the motion of viscous debris flow by using the Bingham Model, which is applicable to laminar viscous flows because in such fluids, the concentration of coarse particles is quite low. Bagnold [9] described it in terms of dilatant flow, and emphasized the discrete force between particles caused by collision. Based on the experiment conducted by Bagnold, Takahashi [10] proposed the idea of the presence of a vertical collision stress between particles, which

supports the coarse particles not sinking into the fluid. Chen [11,12] comprehensively focused on the elastic shear force, the plastic shear force, and the laminar shear force in muddy fluids, implementing a modification of the motion resistance which could be suitable for viscoplastic fluids containing certain coarse particles. O'Brien [13,14] took into account the effects of plasticity, viscosity, collision, and turbulence on resistance forces and identified a general equation which includes the components due to cohesive yield stress, Mohr-Coulomb shear force, viscous shear stress, turbulent shear stress, and dispersive shear force, establishing a debris fluid model that could fully describe the effects due to different particle compositions and different densities.

Furthermore, the mesh-free numerical modeling techniques are being gradually introduced into this field. Dai [15] developed a three-dimensional model to simulate rapid landslide motions across 3D terrains. The artificial viscosities linearly related to the linear and quadratic terms of shear deformation were incorporated into the pressure terms in the momentum equation to dissipate energy for avoiding numerical oscillations. However, Dai [15] did not consider the effect of yield stress or the interaction between solid and liquid phases. Hosseini [16] adopted an innovative treatment similar to the one applied in Eulerian formulations to viscous terms, and hence facilitated the implementation of various inelastic non-Newtonian models. This approach utilized the exact forms of the shear strain rate tensor and its second principal invariant to calculate the shear stress tensor. Rodriguez-Paz [17] introduced a new frictional approach for the boundary conditions and modified constitutive equations in the SPH (Smoothed Particle Hydrodynamics) method to focus on the interaction between fluid particles and boundary conditions. The resulting technique was then applied for the numerical simulation of debris flows and the results were compared with those experimentally obtained and available in literature [18,19], providing good agreements.

In this paper, the SPH method was applied to estimate the movement of a non-homogeneous viscous debris-flow: the fluid under investigation was divided into solid-liquid phases. The solid phase was characterized by particles with larger size than the critical particle size, while the liquid phase was the mixture of water and particles with smaller size than the critical one. The constitutive relation for the liquid phase was characterized by yield stress, laminar viscous force, and turbulent viscous force, while the constitutive relation of the solid phase was related to support force, friction, and collision stresses. It is well known that the magnitude of shear deformation determines which force plays a dominant role in the process of the fluid movement [20]. Therefore, quantifying all the above forces, it could be possible to quantify and estimate the role of each component and further investigate how the shear sharp deformation could be reduced. However, due to the complexity of the materials composing debris flows, even the fluid with same rheological coefficients could generate effects attributable to different flow structures and characteristics. Therefore, the effect on the debris flow process of the initial vertical distribution of the two-phase solid-liquid is also considered in this study.

#### **2. Fundamental Theories and Numerical Modeling**

#### *2.1. The SPH Method*

SPH is a kind of mesh-free method based on a pure Lagrangian description, which has been widely applied in multiple engineering and science fields [21–24]. Compared with the mesh method based on continuum theory, the SPH method avoids the problem of mesh distortion in dealing with the flow issue since there is no connectivity between the particles, since it is developed on a uniform smoothed particle hydrodynamic framework. By adopting this technique, the goal is to provide accurate and stable numerical solutions for integral equations or partial differential equations (PDEs) using a series of arbitrarily distributed particles carrying field variables, such as mass, density, energy, and stress tensors [25].

#### 2.1.1. SPH Interpolation

There are two main steps to transform PDEs equations into the SPH form, called particle approximation and kernel approximation, respectively [26]. The first step consists of representing a function in continuous form as an integral representation by using an interpolation function. In this step, the approximation of the function and its derivatives is based on the evaluation of the smooth kernel function and its derivatives. The second step involves representing the problem domain by using a set of discrete particles within the influence area of the particle at location *x*, and then estimating the field variables for those particles as follows:

$$
\langle f(\mathbf{x}) \rangle = \int\_{\Omega} f(\mathbf{x'}) \mathcal{W}(\mathbf{x} - \mathbf{x'}, h) d\mathbf{x'}, \tag{1}
$$

$$\mathcal{N}\{f(\mathbf{x})\} = \sum\_{j=1}^{N} m\_j \frac{f(\mathbf{x}\_j)}{\rho\_j} \mathcal{W}(\mathbf{x} - \mathbf{x}\_{j'} h) . \tag{2}$$

where *x*' denotes the position of continuum in the influence domain before the discretization; *W* denotes the smoothing function; *h* is a parameter that defines the size of the kernel support, known as the smoothing length; Ω represents the problem space whose radius is taken as several times of *h* according to different smoothing functions; *N* is the total number of neighboring particles; *m* is the mass; and ρ is the density.

Kernel approximation is the technique of approximating the values of both the field function and the derivative of the field function. The kernels used in the SPH method approximate a δ function (the Dirac function). Monaghan [27] suggested that a suitable Kernel approximation must have a compact support in order to ensure zero interactions outside its computational range. The original calculations of Gingold and Monaghan [28] used a Gaussian Kernel. The Gaussian Kernel function is simple to use and has high accuracy. Especially for the case of disordered particle distribution, this technique generates stable and accurate approximation results. However, the Gaussian Kernel function does not have a strict compact support unless the size of the Kernel support approaches the infinity value. Additionally, further various Kernels forms with a compact support (such as spline [29], super-Gaussian [30], polynomial [31], and cosine [32]) were proposed in previous studies but the one of the most popular Kernels more commonly utilized is based on the spline functions [29] as defined by:

$$\mathcal{W}(r,h) = \frac{10}{7\pi h^2} \times \begin{cases} 1 - 1.5q^2 + 0.75q^3 & 0 \le q < 1\\ 0.25\left(2 - q\right)^3 & 1 \le q < 2\\ 0 & 2 \le q \end{cases},\tag{3}$$

where *q* = |*r*|/*h* and *r* is the separation distance between the particles. In this study, Equation (3) has been used to approximate the values of the field under investigation. This Kernel has compact support so that its interactions are exactly zero for *r* > 2*h*. The smoothing distance or so called "Kernel range" *h* determines the degree with which a particle interacts with adjacent particles.

#### 2.1.2. Gradient and Divergence

As a standard procedure, the gradient and divergence operators need to be formulated in a SPH algorithm if the simulation of the Navier-Stokes equations is to be attempted. In this work, the following commonly used forms are employed for gradient of a scalar A and divergence of a vector *A* [33]:

$$\frac{1}{\rho\_a} \nabla\_a A = \sum\_b m\_b \left( \frac{A\_a}{\rho\_a^2} + \frac{A\_b}{\rho\_b^2} \right) \nabla\_a \mathcal{W}\_{ab\prime} \tag{4}$$

$$\frac{1}{\rho\_a} \nabla\_a \cdot \mathbf{A} = \sum\_b m\_b \left( \frac{\mathbf{A}\_a}{\rho\_a^2} + \frac{\mathbf{A}\_b}{\rho\_b^2} \right) \cdot \nabla\_a \mathcal{W}\_{ab\prime} \tag{5}$$

where <sup>∇</sup>*aWab* is the gradient of the Kernel function *<sup>W</sup> x* − *xj*, *h* with respect to *xj*, subscripts <sup>a</sup> and <sup>b</sup> represent the target particles and the particles in the influence domain, respectively, affecting the position of particle *i*. This choice of discretization operators ensures that an exact projection algorithm is produced. To date, there are various options to represent these operators, but only certain specific ones [34,35] have proven to be more convenient in terms of the accuracy and robustness of the method.

#### *2.2. Governing Equations*

The governing equations for transient compressible fluid flow include the conservation of mass and momentum equations. In a Lagrangian framework, these can be written as follows:

$$\frac{1}{\rho} \frac{D\rho}{Dt} + \nabla \cdot \mathbf{v} = 0,\tag{6}$$

$$\frac{D\mathbf{v}}{Dt} = \mathbf{g} + \frac{1}{\rho} \nabla \cdot \mathbf{\boldsymbol{\tau}} - \frac{1}{\rho} \nabla P\_{\prime} \tag{7}$$

where *t* is time, **v** is the particle velocity vector, **g** is the gravitational acceleration, *P* stands for pressure, and *D*/*Dt* refers to the material derivative. The density ρ has been intentionally kept in the equations to be able to enforce the incompressibility of the fluid. Using an appropriate constitutive equation to model the shear stress tensor τ, Equations (6) and (7) can be used to solve both Newtonian and non-Newtonian flows.

The momentum equations include three driving force terms, i.e., body force, forces due to divergence of the stress tensor, and the pressure gradient, and these always have to be considered together with the incompressibility constraint. In a SPH formulation, the above system of governing equations must be solved for each particle at each time-step, and the order with which the force terms are incorporated into the momentum equations can be different from one algorithm to another.

#### 2.2.1. Equation of State

SPH method can be formulated for incompressible and compressible flows. The equation of state, giving the relationship between particle density and fluid pressure, can be written as follows [28]:

$$P = P\_0 \left[ \left( \frac{\rho}{\rho\_0} \right)^7 - 1 \right]. \tag{8}$$

where *P*<sup>0</sup> represents a constant value of pressure, usually expressed in terms of initial pressure and ρ<sup>0</sup> is the reference density.

#### 2.2.2. Viscous Terms

In the context of the SPH method, several forms of viscosity terms were introduced by Lucy [25], Gingold and Monaghan [29], Wood [34], Loewenstein and Mathews [36], and Shao and Lo [37]. As the purpose of this work was to solve non-Newtonian fluids, a new description of viscosity was developed to facilitate the modeling of such flow characteristics. Viscosity of incompressible Newtonian fluids depends only on the second principal invariant of the shear strain rate:

$$\mathbf{D} = \begin{bmatrix} 2\frac{\partial \mathbf{u}}{\partial x} & \frac{\partial \mathbf{u}}{\partial y} + \frac{\partial \mathbf{u}}{\partial x} \\ \frac{\partial \mathbf{u}}{\partial y} + \frac{\partial \mathbf{u}}{\partial x} & 2\frac{\partial \mathbf{u}}{\partial y} \end{bmatrix} \tag{9}$$

*Water* **2019**, *11*, 2314

In solving Equations (8) and (9), the finite difference method should firstly be used to solve the total derivative between two particles, and then decompose the results into *x*, *y* directions, following Lo and Shao [38].

$$
\left(\frac{\partial u\_i}{\partial \mathbf{x}\_j}\right)\_a = \left(\frac{\partial u\_i}{\partial r\_{ab}}\right)\left(\frac{\partial r\_{ab}}{\partial \mathbf{x}\_j}\right) = \frac{(u\_i)\_a - (u\_i)\_b}{r\_{ab}}\frac{\left(\mathbf{x}\_j\right)\_a - \left(\mathbf{x}\_j\right)\_b}{r\_{ab}},\tag{10}
$$

The viscous debris flow studied in this paper is a non-Newtonian fluid whose constitutive equation has the following form:

$$\text{Bingham fluid}: \ \pi = \pi \mathfrak{z} + \mu \mathfrak{y} \mathbf{D},\tag{11}$$

$$\text{dilatant fluid}: \ \tau = \mu\_D(|\mathbf{D}|) \mathbf{D}, \tag{12}$$

$$\text{Viscopelastic fluid}: \ \tau = \frac{1}{2} \left( \frac{\mu\_0}{|\mathbf{D}|} + \mu\_1 + \mu\_2 |\mathbf{D}| \right) \mathbf{D},\tag{13}$$

where τ*<sup>B</sup>* refers to the yield stress; μ*B*, μ<sup>1</sup> are the coefficients of friction and μ*D*, μ<sup>2</sup> are the coefficients of collision stresses. According to the different flows considered, the symbols in Equation (13) can represent different meanings. When the object considered corresponds to solid particles, μ<sup>0</sup> indicates the static support force between the solid particles, μ<sup>1</sup> is the particle friction coefficient, and μ<sup>2</sup> is the particle collision coefficient. When considering a liquid phase slurry, μ<sup>0</sup> is the yield stress, μ<sup>1</sup> is the laminar viscosity coefficient, and μ<sup>2</sup> is the turbulent viscosity coefficient. By comparing Equations (10)–(12), Equations (10) and (11) can be regarded as a special form of Equation (12), because slurry flows consists of water and fine particles (liquid phase) and coarse particles (solid phase). In the Bingham model, since the fluid turbulence is not considered, the turbulent viscosity coefficient is μ<sup>2</sup> = 0. In the expansion model, the inter-particle frictional force is negligible relative to the particles' collision, so the second term in Equation (11) is zero.

By substituting single components of Equations (9) and (10), the second term of right hand side in Equation (7) can be written as follows:

$$\frac{1}{\rho\_d} \nabla\_d \cdot \mathbf{r} = \frac{1}{\rho\_d} \begin{bmatrix} \frac{\partial}{\partial x} \\ \frac{\partial}{\partial y} \end{bmatrix} \cdot \begin{bmatrix} \tau\_{xx} & \tau\_{xy} \\ \tau\_{yx} & \tau\_{yy} \end{bmatrix} = \frac{1}{\rho\_d} \begin{bmatrix} \frac{\partial}{\partial x} \cdot \tau\_{xx} + \frac{\partial}{\partial y} \cdot \tau\_{yx} \\ \frac{\partial}{\partial x} \cdot \tau\_{xy} + \frac{\partial}{\partial y} \cdot \tau\_{yy} \end{bmatrix} \tag{14}$$

By substituting Equation (13) into Equation (14), the viscous term in the *x*, *y* direction can be represented by the following discretization:

$$\frac{du}{dt} = \sum\_{b} m\_{b} \left\{ \frac{1}{\rho\_{a} \rho\_{b}} \cdot \frac{1}{2} \left( \frac{\mu\_{0}}{|\mathbf{D}|} + \mu\_{1} + \mu\_{2} |\mathbf{D}| \right) \Big| 2 \frac{\partial u}{\partial \mathbf{x}} \mathbf{i} + \left( \frac{\partial v}{\partial \mathbf{x}} + \frac{\partial u}{\partial y} \mathbf{j} \right) \Big| \right\} \cdot \left( \frac{\partial \mathcal{W}}{\partial q} \cdot \frac{\mathbf{x}\_{ab}}{hr} \mathbf{i} + \frac{\partial \mathcal{W}}{\partial q} \cdot \frac{y\_{ab}}{hr} \mathbf{j} \right) \tag{15}$$

$$\frac{d\upsilon}{dt} = \sum\_{b} m\_b \left\{ \frac{1}{\rho\_a \rho\_b} \cdot \frac{1}{2} \left( \frac{\mu\_0}{|\mathbf{D}|} + \mu\_1 + \mu\_2 |\mathbf{D}| \right) \left[ \left( \frac{\partial \upsilon}{\partial \mathbf{x}} + \frac{\partial \mu}{\partial y} \mathbf{j} + 2 \frac{\partial \upsilon}{\partial y} \mathbf{j} \right) \right] \cdot \left( \frac{\partial \mathcal{W}}{\partial q} \cdot \frac{\mathbf{x}\_{ab}}{hr} \mathbf{i} + \frac{\partial \mathcal{W}}{\partial q} \cdot \frac{y\_{ab}}{hr} \mathbf{j} \right) \tag{16}$$

It can be seen from Equation (15) that when deformation |**D**| is relatively small, yield stress (particle support force) has a greater impact on the fluid's acceleration, but there should be an upper limit to this effect in order to prevent excessive acceleration which could cause the local instability of the fluid investigated. In existing studies [16], a lower limit is usually set for |**D**|. When |**D**| is less than this lower limit, the relationship between stress and |**D**| satisfies linearity:

$$\begin{aligned} |\mathbf{D}| \le \frac{\mu\_0}{\sigma} \to \tau = \sigma \mathbf{D} \\ |\mathbf{D}| > \frac{\mu\_0}{\sigma} \to \tau = \frac{1}{2} \left( \frac{\mu\_0}{|\mathbf{D}|} + \mu\_1 + \mu\_2 |\mathbf{D}| \right) \mathbf{D} \end{aligned} \tag{17}$$

where σ is the limiting factor.

In [16], Hosseini considers that the viscosity of "solid zone" fluid is much greater than that of the main fluid (100 times). In this study, although the turbulent stress term <sup>1</sup> <sup>2</sup>μ2|**D**|**D** is introduced, it can be ignored in the region of low velocity. By calibrating this value against experimental results obtained for this study, 200μ<sup>1</sup> was the limited factor selected.

#### **3. Experimental Setup and Boundary Conditions**

Figure 1 shows the facility where the experiments for the numerical validation of the method previously described in Section 2 were conducted. This facility is located in the Debris Flow Observation and Research Station at Jiangjia Gully, the largest field research center in China—also known as "the debris flow museum". For this study, water is stored upstream to the material prepared under multiple configurations and water is always released by opening a steel gate at the velocity of 3 m/s. The flume can then be divided into two parts: (i) a steep upstream reach (6.0 m long and the chute slope can vary from 15◦ to 40◦ (always set up at 25◦ for these experiments); and (ii) a flat-bottomed downstream section (3.0 m long). The slope of the flume's bed can be manually adjusted.

Details regarding the experimental procedure conducted can be found in [6]. For this study, three slurries were testes with densities ρ = 1400 kg/m3; ρ = 1500 kg/m<sup>3</sup> and ρ = 1600 kg/m3. Different layers patterns were selected according to the different configurations displayed in Figure 2. By adding water to the flume, when the water level reached the height of the mixing fluid, the front-end steel gate of the mixing area was released at a speed of 3 m/s. In order to maintain the driving force of the mixtures, the water level behind the mixtures was kept at *h* = 0.2 m during the experiment.

**Figure 1.** Experimental equipment for debris flow simulation [6].

The experimental configurations chosen to examine the influence of vertical grading patterns on the debris flow formation and propagation are displayed in Figure 2.

**Figure 2.** Schematic representation (experimental and numerical) of configurations applied for this study [6].

There are three configurations: Figure 2a shows the configuration a, with coarse particles as a top layer and fine sediment positioned below them; Figure 2b shows the configuration b, a realistic configuration with the original sediments collected from the Jiangjiagou Valley, fully mixed and distributed along the entire measurement area; Figure 2c shows the configuration c, with coarse particles distributed at the bottom and fine grains on the top of them. For this study, three types of liquid slurries were considered: (i) ρ = 1400 kg/m3, (ii) ρ = 1500 kg/m3, and (iii) ρ = 1600 kg/m3. The slurry rheology coefficient was measured by the MCR301 advanced rotary rheometer manufactured by Anton Paar, Austria [39]. Values of viscosity μ for the fluids simulated are displayed in Table 1.


**Table 1.** Experimental conditions (viscosity values and vertical grading patterns) adopted for this study.

The particle size distribution curves are shown in Figure 3, Ψ stands for the percentage of accumulated mass of particles and *d* denotes particle size.

**Figure 3.** Particle size distribution curves.

#### **4. Simulation and Analysis**

By simulating the movement processes of different viscous debris flows, a series of experiments has been completed where free surface heights, fluid velocities, pressures, and shear deformations associated with the movement of the fluid were measured. The numerical simulation was carried out to generate the experimental conditions previously described (Figure 4). The numerical simulation replicated *N* = 113,054 particles, particle spacing *dp* = 0.0025 m, solid particle density ρ*<sup>s</sup>* = 2200 kg/m3, thickness of solid phase *hs* = 0.1 m and thickness of liquid phase *hl* = 0.1 m for configuration a and c in Figure 2. Similar to the tests conducted on the experimental facility, three different viscosity coefficients for the liquid phases (as shown in Table 1) were selected in the numerical simulation. The inflow conditions were the same as those applied experimentally, and the water level as driving force of the debris flow was kept at 0.2 m for each entire simulation.

**Figure 4.** Initial state of debris flow numerical simulation (mixed configuration). The water level behind the mixtures was kept at *h* = 0.2 m.

Figure 5 represents the free surface values recorded at different times for each of the tests conducted in Table 1, after the debris flow initiation generated by the release of water upstream. The *x*-axis represents the length *L* of the debris flow, and the *y*-axis represents the height *H* of the debris flow. After 0.7 s, the distance reached by the debris flow (considering all the configurations) is within the range 2.48–2.66 m, and the maximum velocity range is 4.65–5.12 m/s. Authors noticed that when the debris flow has similar viscous properties but for the vertical distribution of different particles, the differences in velocities are not so significant and can considered almost negligible in most of the cases. However, the maximum velocity recorded for configuration b (Shown in Figure 2b correspond to tests 2, 5 and 8 in Figure 5) is similar to the one recorded for configuration a (Shown in Figure 2a and corresponding to tests 1, 4, and 7 in Figure 5), while configuration c (Shown in Figure 2c and corresponding to tests 3, 6, and 9 in Figure 5) was characterized by higher values of velocities and elevations measured. By comparing the shapes of head under different vertical distributions, it was found that for the tests conducted in Table 1, free surface values measured for configuration b fluctuate more than in configuration a and c, demonstrating that this scenario is typical of intermittent debris flows.

**Figure 5.** Debris flow free surface for tests 1, 2, 3, 4, 5, 6, 7, 8, and 9 (Table 1) and experimental results used for validation (dots), respectively.

#### *Characterization of Intermittent Debris Flow*

In order to study the causes of this phenomenon, the characteristics of the fluid movement processes associated to configuration b were analyzed. Figure 6 shows the evolution of the solid-liquid phase at different locations simulated numerically. Figure 6a shows that in the horizontal region, most of the particles still retain under the laminar form. When moving into the upstream part of the sloped section, the fluid height decreases and the liquid phase group is stretched, as shown in Figure 6b. Then, due to the slope, velocity increases while the fluid height decreases, and different layers of liquid and solid particles will appear almost as parallel mixing within the entire width of the debris flow, as shown in Figure 6c. At this stage, the altering layers interact changing continuously positions

demonstrating that mixing processes are taking place and when the mixing is finally completed, the fluctuation amplitude reduce becoming more stable, while the influenced range of the fluctuations can spread over a longer length, as shown in Figure 6d. Figure 6e shows the effect of the liquid phase on the height of the debris flow. It is clear that when liquid particles accumulate due to the mixing phenomena (highlighted as circles in Figure 6e), there is a correspondent decrease of the height of the debris flow (pointed out using arrows in Figure 6e). This inverse relationship is very interesting especially because it demonstrates how the gathering and accumulation of liquid particles tends to appear towards the bottom side of the debris flow layer.

(**e**)

**Figure 6.** Analysis of the debris flow behaviors at different locations. Solid particles and liquid particles are represented by red and blue dots, respectively.

On this basis, the relationship between the moisture content φ (the amount of liquid particles divided by the amount of solid particles, *Nl*/*Ns*), the kinetic energy of particles *Ek* and the height of the free surface *H* (related to the potential energy of particles *Ep*) were calculated for the tests conducted for configuration b. As shown in Figure 7a, the height of the free surface *H* decreases as the debris flow develops. There is a noticeable correlation between the fluctuation of the free surface associated with the fluctuation of the moisture content. In the regions of *L* = 1.00–1.38 m and *L* = 1.82−2.04 m, the height of free surface decreases linearly, and in these two regions, the water content remains in the range 0.2–0.65. The points that obviously exceed this threshold are *L* = 0.90, *L* = 1.52, *L* = 1.6, *L* = 1.74, and *L* = 1.80−1.84, and the height of free surface is different from that of linear decline in these areas or vicinity. When the moisture content is within the range 0.20–0.65, the free surface of the debris flow is characterized by a linear change, but when the moisture content exceeds this range, it generates an impact on the free surface.

**Figure 7.** Relationships between the moisture content, the kinetic energy, and the height of the free surface for configuration b.

From Figure 7b, it can be seen that there is a more obvious negative correlation between the particles kinetic energy *Ek* and the moisture content φ. To almost every peak of the kinetic energy *Ek* (highlighted as green circles in Figure 7b) calculated corresponds a peak of the moisture content φ (highlighted as blue circles in Figure 7b), which indicates that kinetic energy *Ek* and moisture content φ interact directly. However, this effect can only be assigned to small-scale portions of the particle kinetic energy fluctuations. Looking at Figure 7b, at the location of *L* = 1.74 m, the moisture content value corresponds to φ = 0.7619 and it is the maximum value measured in this region, and the corresponding kinetic energy *Ek* records its minimum value. But because of the large kinetic energy of the particles recorded in this region, the corresponding kinetic energy *Ek* = 5.4848 J is still higher than that recorded at the position of *L* = 1.64 m in the adjacent one *Ek* = 0.30263 J.

Figure 8 shows the effect of moisture content φ on the kinetic energy *Ek* and the height of the free surface *H* on a large scale, for configurations a and c. For both configurations, where the solid phase is located at the top and the bottom, the free surface is greatly affected by the magnitude of the moisture content φ, while the kinetic energy *Ek* is greatly affected by the derivative of moisture content along the length of the slope *<sup>d</sup>*<sup>ϕ</sup> *dL* . However, the fluctuation of the moisture content Δφ along the length *L*, especially for configuration a where the solid phase is displayed at the top of the debris flow, is relatively small. So the kinetic energy *Ek* and potential energy *Ep* curves show relatively large-scale area fluctuations and linear characteristics in comparison to the mixed distribution fluid conditions typical of configuration b.

**Figure 8.** Relationships between the moisture content φ, the kinetic energy *Ek*, and the height of the free surface *H* for configuration a and c.

Finally, the energy conversion curves of fluids with different viscous coefficients were inspected and confronted, as shown in Figure 9. It was found that the gravitational potential energy (*Ep* = *mgH*) and the total energy (*E*<sup>0</sup> = *Ek* + *Ep*) of fluids decreases at a similar rate. The difference between three

fluids is mainly reflected on kinetic energies. When comparing the set of fluids with the smallest density (ρ = 1400 kg/m3, μ<sup>0</sup> = 0.00004, μ<sup>1</sup> = 0.0048, and μ<sup>2</sup> = 0.0197), results shows that velocity values increase from *t* = 0.0 s up to *t* = 0.6 s, reaching almost the highest values, and then the kinetic energy of the three fluids tends to be equal. As the time progresses, the same order appears again in the kinetic energy magnitude arrangement, which is *Ek,*<sup>ρ</sup> = 1400 kg/m<sup>3</sup> > *Ek,*<sup>ρ</sup> = 1500 kg/m3 > *Ek,*<sup>ρ</sup> = 1600 kg/m3. This phenomenon is also due to the stronger fluctuation of the less dense fluids and these effects caused by different viscous fluids on debris flow array and collision, and friction forces on debris flow movement, will require further investigation in the future.

**Figure 9.** Energy evolution of different bulk density debris flows.

#### **5. Summary and Conclusions**

Debris flows are a natural phenomenon causing a lot of economical and human losses worldwide. Due to their nature, debris flows travel long distances at high speeds, and the time-space evolution of the relationship between soil and water content strongly affects the propagation stage. Thus, a quantitative modeling of this phenomenon is crucial to design strategies to be adopted to reduce the negative impacts. This paper contributes to this topic through the use of a SPH model investigating multiple combinations of fine and coarse particles with water content.

The SPH model was used to simulate tests conducted in the experimental facility that is located in the Debris Flow Observation and Research Station at Jiangjia Gully, the largest field research center in China—also known as "the debris flow museum". As previously demonstrated, the SPH model is capable of properly reproducing the main characteristics of debris flows (propagation height and velocity, and more importantly to correctly simulate the time-space evolution of solid and liquid particles during the whole process from initiation to propagation over an impervious/permeable bottom boundary).

Based on the theory of solid-liquid two phase flows, the viscous term in the SPH model was modified to make it suitable for nonhomogeneous viscous debris flows. It was found that the denser the fluid is, the greater are the yield stress and the turbulent viscous coefficient. However for laminar viscous coefficients, the fluid with the density of ρ = 1500 kg/m<sup>3</sup> has the largest values. The results obtained can be summarized as follows:


The authors can also confirm that there are still some uncertainties within the results analyzed that could be reduced by the use of either novel physically-based entrainment laws or fully 3D mathematical approaches, which could surely more accurately take into consideration the variation of pore water pressures inside the propagating mass. Therefore, future research will target the development of 3D mathematical models to refine the findings and provide an even better understanding of this very complex natural phenomenon.

**Author Contributions:** All the authors jointly contributed to this research. A.S. was responsible for the proposition and design of the experiments, analysis of the results and conclusions of the paper; S.W. completed the numerical simulation, S.W. and M.R. analysed the experimental datasets and wrote the paper. M.W. and J.Q. performed the experiments.

**Funding:** The research reported in this manuscript is funded by the National key research and development plan (Grant No. 2018YFC1406404), the Interdiscipline Research Funds of Beijing Normal University and the Natural Science Foundation of China (Grant No. 11372048).

**Conflicts of Interest:** The authors declare no conflict of interest.

#### **References**


© 2019 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (http://creativecommons.org/licenses/by/4.0/).

### *Article* **Comparative Study on Violent Sloshing with Water Jet Flows by Using the ISPH Method**

#### **Hua Jiang 1,2, Yi You 1, Zhenhong Hu 1,\*, Xing Zheng <sup>1</sup> and Qingwei Ma 1,3**


Received: 13 October 2019; Accepted: 3 December 2019; Published: 9 December 2019

**Abstract:** The smoothed particle hydrodynamics (SPH) method has been playing a more and more important role in violent flow simulations since it is easy to deal with the large deformation and breaking flows from its Lagrangian particle characteristics. In this paper, the incompressible SPH (ISPH) method was used to simulate the liquid sloshing in a 2D tank with water jet flows. The study compares the liquid sloshing under different water jet conditions to analyze the effects of the excitation frequency and the water jet on impact pressure. The results demonstrate that the water jet flows can significantly affect the impact pressures on the wall caused by violent sloshing. The main purpose of the paper is to test the ISPH ability for this study and some useful regulars that are obtained from different numerical cases and study the effect of their practical importance.

**Keywords:** ISPH; liquid sloshing; water jet flow; impact pressure; excitation frequency

#### **1. Introduction**

The phenomenon of violent sloshing appears widely in the field of Naval Architecture and Ocean Engineering, especially in liquid cargo carriers, such as the LNG (liquefied natural gas) carriers. During this process, the motion of liquid in partially filled tanks may cause large global and local loads on the tank walls when the frequency of sloshing is close to the natural frequency of the liquid tank. This consequence would be very serious in engineering practice, which may cause damage to the hull structure and even affect the stability of the carrier [1]. Therefore, sloshing is an old topic, but it still needs to be studied in depth.

At present, lots of studies on the violent sloshing flows have been carried out by the linear and nonlinear potential flow theory and the scaled model experiment. Faltinsen [2–4] used the incompressible potential flow theory to simulate liquid sloshing and obtained the formulas that have been widely used in the field of sloshing simulation. However, the method can be used to study sloshing tanks with relatively simple geometry and internal structure. In addition, Akyildiz et al. [5] investigated the pressure distribution on a rectangular tank during the process of sloshing by an experimental method. Sames et al. [6] studied sloshing in a rectangular tank with a baffle, and a cylindrical tank was also considered. Indeed, the experimental method can be applied to study the sloshing in the tank with more complex shapes, but it also requires high expenses for the site and facility. Hence, the numerical method has been getting more important in the simulation of liquid sloshing in recent years. The conventional numerical methods are carried out by using Euler grids. Wu et al. [7] simulated the sloshing waves in a 3D tank based on the finite element method (FEM). In the conventional grid-based methods, in order to track the moving free surface, some additional

techniques, such as the Volume-of-Fluids, are used in the methods. The VOF uses the volume fraction of fluid in gird to define the free surface. However, the problems of numerical diffusion become serious when the surface cell becomes extremely complicated, such as in a liquid sloshing, which can easily fail to simulate because of the large deformation of grids. Recently, a kind of mesh-less method named smoothed particle hydrodynamics (SPH) has attracted quite a few researchers' attention [8]. It does not depend on any grids, and the computation is purely based on a group of discrete points that can move freely. So, it can capture the free surface flow conveniently, which is more suitable for treating the problems of large deformation of free-surfaces. Delorme and Colagrossi et al. [9] investigated impact pressure in the case of shallow water sloshing by the SPH method, compared the results with experimental ones, and then discussed the influence of viscosity and density re-initialization on the SPH results. Gotoh and Khayyer [10] simulated the violent sloshing flows using the incompressible SPH (ISPH) method and presented two schemes to enhance the accuracy of the simulation of impact pressures. Zheng and You [11] compared the effect of different baffle configurations on mitigating sloshing by the ISPH method. A great deal of research [12,13] shows that the ISPH method can improve the accuracy and stability of the calculation pressure, and the pressure field is smoother.

As a matter of fact, the marine environment is always very complex when sloshing happens. If an oil fire occurs, it would be a big disaster and hard to control, and the water jet flow outside would get into the tank to put out the oil tank fire, which would influence the impact effect of sloshing. Such a consequence may cause more serious damage to the hull structure, which would be very serious in engineering practice. With regard as the problem of jet flows, Hatton et al. [14,15] studied the trajectories of large water jets that are used in the design of fire-fighting systems, particularly those used in offshore situations, and evaluated the effects of flow-rate, pressure, and nozzle size during the process of the system design. Fischer et al. [16] used three different CFD codes, namely, the CHYMES multiphase flow model, the FEAT finite element code, and the Harwell-FLOw3D finite volume code, to simulate the problem of a laminar jet of fluid injected into a tank of fluid at rest and make a detailed comparison. Aristodemo et al. [17] studied the plane jets propagating into still fluid tanks and current flows by using the WCSPH method. Andreopoulos et al. [18] carried out an experiment on the flow generated by a plane with a buoyant jet discharging vertically into shallow water.

In this paper, the liquid sloshing with a water jet flow from the top of the tank will be studied by using the incompressible SPH (ISPH) method. Through the comparison of different situations, the sloshing effects and characteristics of the impact pressure are studied. The aim of this study is to summarize the influence of the water jet flow on sloshing, so as to give a reference for practical engineering.

#### **2. ISPH Methodology**

#### *2.1. Governing Equations*

The SPH model is based on the semi-Lagrangian form of the continuity equation and the momentum equation. In the ISPH method, the density of fluid is considered to be a constant, and thus, the governing equations are written as follows:

$$
\nabla \cdot \mathfrak{u} = 0,
$$

$$\frac{Du}{Dt} = -\frac{1}{\rho}\nabla P + \mathbf{g} + \nu\_0 \nabla^2 u\_\prime \tag{2}$$

where ρ is the density of fluid; *u* is the velocity of particle; *t* is the time; *P* is the particle pressure; *g* is the gravitational acceleration; ν<sup>0</sup> is the kinematic viscosity; and ∇ is Hamilton operator, which is a vector operator.

#### *2.2. Particle Approximation*

The computational domain of the SPH method is composed of a group of discrete particles, and each particle is given corresponding physical information, such as density, volume, mass, velocity, and pressure. The physical information of each particle can be approximately obtained by the information carried by the surrounding particles, which is shown as follows.

$$f(r\_i) = \sum\_{j=1}^{N} \frac{m\_j}{\rho\_j} f\left(r\_j\right) W\left(r\_{ij}\right) \tag{3}$$

where *f(r)* represents the physical information of particles, *m* is the mass of the particle, *i* and *j* are the center particle and neighbor particle, respectively. *N* is the number of neighbor particles. *W(rij)* is the kernel function, which can reflect the different effects between different particles. In this paper, the cubic B-spline kernel proposed by Monaghan et al. [19] is used as follows:

$$\mathcal{W}\left(r\_{ij},\ h\right) = \alpha\_d \begin{cases} \frac{2}{3} - q^2 + \frac{1}{2}q^3, & 0 \le q < 1\\ \frac{1}{6} \left(2 - q\right)^3, & 1 \le q < 2\\ & 0, \ 2 \le q \end{cases} \tag{4}$$

where *h* is the kernel smoothing length, *rij* is the distance between the *i* and *j* particle, α*<sup>d</sup>* is a constant, and when the case is 2D, its value is <sup>15</sup> <sup>7</sup>π*h*<sup>2</sup> and *<sup>q</sup>* <sup>=</sup> *<sup>r</sup>ij h* .

So the derivatives of *f(r)* can be represented as:

$$\nabla f(\mathbf{r}\_i) = \sum\_{j=1}^{N} \frac{m\_j}{\rho\_j} f\left(\mathbf{r}\_j\right) \nabla\_i \mathcal{W}\left(\mathbf{r}\_{ij}\right),\tag{5}$$

where ∇*<sup>i</sup>* is the gradient, which is taken with respect to the particle *i*.

#### *2.3. Poisson Pressure Equation*

In the ISPH, a two-step projection method is used to solve the velocity and pressure field from the continuity equation and momentum equation [20]. The first step is the prediction of velocity without considering the pressure term. The second step is the correction step in which the pressure term is added through the pressure Poisson equation (PPE), then the PPE is obtained as follows:

$$
\nabla^2 P^{t + \Delta t} = \frac{\rho \nabla \cdot \mathbf{u}^\*}{\Delta t} \tag{6}
$$

where *u\** is the intermediate particle velocity at the first step.

Similarly, Shao and Lo [20] proposed a projection-based incompressible approach by imposing the density invariance on each particle, leading to the following PPE equation:

$$\nabla \cdot (\frac{1}{\rho^\*} \nabla P^{t+\Delta t}) = \frac{\rho\_0 - \rho^\*}{\rho\_0 \Delta t^2} \tag{7}$$

where ρ<sup>∗</sup> is the density at the intermediate time step, ρ<sup>0</sup> is the initial fluid density, and the combined PPE incorporates both the velocity-divergence-free condition and the zero-density-variation condition, which is obtained as:

$$
\nabla^2 P^{t+\Delta t} = a \frac{\rho\_0 - \rho^\*}{\Delta t^2} + (1 - a) \frac{\rho\_0 \nabla \cdot \mathbf{u}^\*}{\Delta t} \tag{8}
$$

where α is a blending coefficient. If α is equal to 1, in Equation (8), the source term of PPE adopts the density variable effect, which may lead to substantial pressure noises and particle randomness caused by larger density changes. If α is equal to 0.0, the source term of PPE adopts the velocity divergence effects, which is smoother for source term distribution, but it will cause the pattern distribution of some particles. In order to make the source term more reliable, α is equal to 0.01 in this paper, according to lots of computational experience. And more advanced PPEs with the error-compensating source term (ECS) can follow the references [21,22].

#### *2.4. Calculation of Spatial Derivatives*

According to Equation (5), the spatial derivatives of pressure and velocity can be calculated as follows:

$$\nabla P\_i = \rho\_i \sum\_{j=1}^{N} m\_j (\frac{P\_j}{\rho\_j^2} + \frac{P\_i}{\rho\_i^2}) \nabla\_i \mathcal{W} \ (r\_{ij}, h), \tag{9}$$

$$\nabla \mathbf{u}\_i = -\frac{1}{\rho\_i} \sum\_{j=1}^{N} m\_j (\mathbf{u}\_i - \mathbf{u}\_j) \nabla\_i \mathcal{W} \ (\mathbf{r}\_{\mathbf{i}j}, h). \tag{10}$$

Therefore, in the ISPH method, the viscous term adopts the following form:

$$\nabla(\boldsymbol{\nu}\_{i}\nabla\cdot\boldsymbol{\mu}\_{i}) = \sum\_{j=1}^{N} 4m\_{j}(\frac{\boldsymbol{\nu}\_{i}+\boldsymbol{\nu}\_{j}}{\rho\_{j}+\rho\_{i}} + \frac{\boldsymbol{\mu}\_{ij}\boldsymbol{r}\_{ij}}{\boldsymbol{r}\_{ij}^{2}+\eta^{2}})\nabla\_{i}\mathcal{W}\left(\boldsymbol{r}\_{ij},h\right) \tag{11}$$

where η is a small parameter to avoid the singularity, and in this paper, the value is chosen to be 0.1 *h*. So, the PPE is discretized by combining the SPH gradient and divergence rules to obtain:

$$\nabla(\frac{1}{\rho^\*}\nabla P) = \sum\_{j=1}^N m\_j(\frac{8}{\left(\rho\_j + \rho\_i\right)^2} \frac{P\_{ij}}{r\_{ij}^{-2} + \eta^2}) \nabla\_i \mathcal{W}\left(r\_{ij}, h\right). \tag{12}$$

The treatment of free surface and solid boundary conditions follows the study of Zheng et al. [23].

#### *2.5. Inlet Boundary Treatment*

With regards to the case of a general closed boundary, sufficient particles are given at the initial time, and no particles are added during the middle steps. However, the model of sloshing with a water jet flow needs to add particles constantly at the top boundary. Hence, the condition of the top boundary should be treated as an inlet boundary to add the particles of the water jet. However, if the particles are added directly, there will be only a row of particles at the beginning, which would lead to big errors in the particle approximation.

In this paper, three rows of virtual particles are arranged at the top boundary, as shown in Figure 1a, which have the same velocity and physical information as the initial water jet particles, but they do not participate in the solving of PPE, they are just used in the particle approximation. The velocities of virtual particles remain constant, and when they have moved a distance of a particle size *dx,* as shown in Figure 1b, they will return to the initial position automatically, as shown in Figure 1c. So, the distance between virtual particles and the water particles is just enough to add a row of new water particles. Finally, the new row of water particles will be added, as shown in Figure 1d.

**Figure 1.** Different steps of the inlet boundary treatment: (**a**) Initial distribution of particles; (**b**) Particles has moved a distance of a particle size; (**c**) The virtual particles are return to the initial positions; (**d**) The new water particles are added.

#### **3. Numerical Results and Validation**

In this section, the numerical results of the sloshing tank and water jet calculated by the ISPH method are compared with the experimental data to validate the accuracy of the ISPH model.

#### *3.1. Sloshing Tank Simulation and Validation*

Figure 2 illustrates the rectangular sloshing tank that was used in the experiment of Liao and Hu [24], and the schematic view gives the geometrical dimensions of the sloshing tank; the length of the tank is 1.2 m, the height is 0.6 m, the initial depth of the liquid inside is 0.12 m, and the rotation center is located at the geometric center of the tank. M1 and M2 are two pressure sensors, which are used to record the pressure of the two points. In the experiment, the sloshing tank is imposed a rolling motion, and the amplitude and period of excitation are set as 10◦and 1.85 s, respectively. The distribution of the particles in the ISPH computational model is uniform, with a particle spacing *dx* = 0.002 m. A constant time step of *dt* = 0.0003 s is used. The pressure data is extracted for 8 s after about twice the sloshing period *T,* which is relatively stable.

Figure 3 shows the snapshots of the sloshing motion of the ISPH model with the calculated pressure fields. Meanwhile, the corresponding experimental photographs are considered for the comparisons of the motion patterns. According to the comparisons, the ISPH model can get a good agreement on surface profile with the experimental results, and most of the features of the violent sloshing process have been captured by the ISPH model in a satisfactory manner. In addition, the ISPH method also provides a reliable regular distribution of pressure at the impact regions, and the pressure fields show a very stable pattern with little pressure noise.

**Figure 2.** Schematic view of the sloshing tank following Liao and Hu [24].

**Figure 3.** Sloshing snapshots computed by the ISPH (left) compared with the experimental photos (middle), and overlaid free surfaces (right): (**a**) 0.0*T*; (**b**) *T*/6; (**c**) *T*/3; (**d**) *T*/2.

Figure 4 presents the time histories of calculated pressure at points M1 and M2 by the ISPH method, which are compared with the experimental data. From the contrast, the pressure traces with both first and second peaks calculated by the ISPH method have made good accuracy with experimental results. So the comparisons between the ISPH results and experimental results demonstrate that the ISPH method can get accurate results for the maximum peak values and the phases of the pressure time histories.

**Figure 4.** Comparisons of pressure time history between ISPH and experimental data of Liao and Hu [24] at: (**a**) M1; (**b**) M2.

#### *3.2. Convergence Analysis of ISPH Model*

In order to validate the convergence of the ISPH sloshing model, three different particle sizes have been tested, which are 0.002, 0.003, and 0.004 m, corresponding to the number of particles are 49,000, 16,000, and 9000, respectively. And the pressure at M1 and M2 calculated in three different particle sizes are compared with the experimental data, and the results are shown in Figure 5. It shows that the amplitude of non-physical pressure oscillations decreases when particle size is from 0.004 m to 0.002 m. Moreover, the enlarged parts manifest that the pressure traces are more smoothed and more closed to the experimental results. Furthermore, the quantitative comparisons in three different particle sizes are carried on M1 and M2, which obtained the values of the mean error in Equation (13), the results are shown in Table 1, more details can refer to Zheng et al. [23]. From the comparison of pressure results, the errors decrease with the particle number increasing, and the results demonstrate that the ISPH method has very good convergence and stability.

$$E\_d = \frac{1}{N} \sum\_{i=1}^{N} \left| \frac{\widehat{P(t)} - P(t)}{P(t)} \right| \tag{13}$$

where *Ea* is the defined mean error; *P* -(*t*) is the numerical results of the pressure; *P*(*t*) is the experimental pressure; *N* is the number of sampling points.

**Figure 5.** Comparisons of pressure time histories between incompressible smoothed particle hydrodynamics (ISPH) and experimental data of Liao and Hu [24] at (**a**) M1 and (**b**) M2.

**Table 1.** Numerical error of ISPH in three particle sizes.


Figure 6 gives the convergence rate of the pressure results calculated by this ISPH method. It is shown that the numerical error *Ea* is closer to the second-order accuracy for the results of pressure calculated at M1 and M2. The ISPH method provides good performance for liquid tank sloshing; thus, it will be used for the study of the sloshing tank applications.

**Figure 6.** Convergence and error analysis of pressure calculated by the ISPH method.

#### *3.3. Numerical Validation of the Water Jet Model*

Figure 7 shows the schematic view of the water jet test case, which was used in the experiment of Kvicinsky et al. [25]. The diameter of the jet inlet is *D* = 0.03 m, and it is located at *H* = 0.1 m above a flat plate. The initial velocity of the water jet was *v* = 19.81 m/s. In the ISPH, a particle size *dx* = 0.0015 m was used, the time step *dt* = 0.00001 s, and the physical time of simulation was *T* = 0.03 s. The pressures at the plate were recorded.

**Figure 7.** Schematic view of the water jet, following Kvicinsky et al. [25].

In the experiment, the pressure coefficient *Cp* on the flat plate was used as follows:

$$C\_p = \frac{2P}{\rho v^2} \tag{14}$$

where *P* is the pressure at the flat plate, ρ is water density, and *v* is the initial velocity of the water jet. Figure 8 shows the snapshots of the numerical results. In order to get the stable pressure data, the pressures were extracted when the physical time is greater than 0.015 s. Figure 9 shows the comparison of pressure coefficient *Cp* among the ISPH results, the experimental data, and the CFD results, where *x* is the coordinate value on the flat plate. Figure 9 shows that the ISPH results have good agreement with both the experimental data and CFD results, which demonstrates that the ISPH method can have high accuracy on the impact pressure caused by the water jet flow.

**Figure 8.** Pressure contour snapshots of the water jet case by ISPH: (**a**) *T*/10; (**b**) *T*/4; (**c**) *T*/3; (**d**) *T*.

**Figure 9.** Comparisons of pressure coefficient *Cp* among the ISPH results, the CFD results, and the experimental data of Kvicinsky et al. [25].

#### *3.4. Validation of Injected Water Jet Flow Model*

Figure 10 gives the geometric dimensions of a rectangular tank. The height of the tank is *H* = 0.5 m, its width is *B* = 0.3 m, and the tank is full of water. *Di* = 0.08 m, which is the inflow jet diameter in the top center of the tank, and *D0* = 0.012 m, which is the outflow diameter. The water jet is injected from *Di* with the velocity *V* = −0.45 *j* m/s, and its direction is straight down.

**Figure 10.** Sketch of the jet injected at the top of a water tank.

Figure 11 shows the comparisons of the vertical component of the velocity field *w* between ISPH and the results of Aristodemo et al. [17] at three significant time instants. It can be seen from Figure 11 that the still water is accelerated progressively with the entry of the water jet flow, and the highest negative velocities are symmetrical by the jet centerline at the initial time. However, these negative velocities become non-symmetrical as time goes on. According to the results of the vertical velocity component by the ISPH method, it is in good agreement with the results of Aristodemo et al. [17].

To quantify the accuracy of the ISPH results, Figure 12 gives the comparisons of relative depth *h*/*H* of jet flow between the results of Aristodemo et al. and Fletcher et al. [16,17], where *h* is the depth of jet flow penetration. The results of ISPH have a high degree of coincidence with the ones of other numerical models. The ISPH method has good potentials for the injected water jet flow.

**Figure 11.** *Cont.*

**Figure 11.** Comparisons of the vertical component of the velocity field between the results of the ISPH method (left) and the results of Aristodemo et al. [26] (right) at *t(V*/*Di)* = (**a**) 2.8; (**b**) 6.2; (**c**) 13.4.

**Figure 12.** Comparison of the relative depth of jet flow penetration between ISPH and other numerical models.

#### **4. Results and Analyses**

The main scientific problem discussed in this paper is the effect of the water jet on sloshing. The calculation formula of the natural frequency of a liquid tank is shown as follows:

$$
\Omega\_0 = \sqrt{(\text{g\pi}/L)\tanh(\pi d/L)}\tag{15}
$$

where *g* is the gravitational acceleration, *L* is the length of the tank, and *d* is the depth of water.

#### *4.1. Sloshing Behaviors with Water Jet*

In this section, the liquid sloshing in a rectangular tank with a water jet from the top of the tank is considered. The configuration of the sloshing tank is shown in Figure 13, and the geometrical dimensions of the sloshing tank are the same as the tank in Figure 2, but one difference is that there is water flow jetted from the center on the top of the tank. The tank experienced a rolling motion with an amplitude of 8◦ and an excitation period of 1.85 s. Four sensors, P1–P4, were placed, as shown in the figure, to monitor the pressure distributions. The initial water depth is 0.12 m, and the tank is rotated at the geometric center. The water jet enters the liquid tank through the center on top of the tank, and the initial velocity was imposed to 0.3 m/s.

**Figure 13.** Schematic view of the sloshing tank with the water jet from the center on top.

The computed particle snapshots with pressure contours and streamlines for this case are shown in Figure 14. It shows the violent sloshing flow and water jet flow with strong interaction for a sloshing period *T*. The smooth and noise-free stable pressure fields indicate the robustness of the presented ISPH method, which is a good way to analyse the impact pressures. It also can be seen from the streamlines that the fluid is entrapped generating some recirculating counter-rotating cells in the violent sloshing process.

Then, the pressure histories at four different locations at the bottom are recorded and compared. Figure 15 gives the comparison of four pressure histories when the velocity is 0.3 m/s, it shows that the trends of four pressure histories are generally similar, the peak value of impact pressure gradually increases with the water jets into the tank, and the two-peak pressure patterns of four histories are remarkable. By way of contrast, the maximum pressure appears at the location of P1, meanwhile the change range of impact pressure is also the largest. It can get higher pressure when the monitor location is closer to the side wall, the sloshing effect is also more violent in this area. Therefore, in the next cases, the impact pressure at P1 will be studied and analyzed.

The liquid sloshing in a rectangular tank with a water jet of various velocities from the top of the tank is considered. The configuration of the sloshing tank is shown in Figure 13, and different initial velocities of the water jet flow are used in this case, which are 0.2, 0.3, 0.4, and 0.5 m/s, respectively. Figure 16 portrays the snapshots of the sloshing process with pressure fields at *t* = 0.296 s with different initial velocities. It can be seen that when the water jet flow gets into the free surface, the strong impact generates a large free surface deformation with two waves running out along the water jet flow, creating an open-air cavity. Then, the air cavity gradually closes, and the two free surfaces form a short-term bump. From the comparison, the open-air cavity becomes larger with the initial velocity increasing.

**Figure 14.** Particle snapshots with pressure contours (left) and streamlines (right) of sloshing flow interactions with water jet flow: (**a**) *T*/4; (**b**) *T*/2; (**c**) 3*T*/4; (**d**) *T*.

**Figure 15.** Comparison of pressure histories at four different locations at the bottom.

**Figure 16.** Comparison of sloshing patterns with different initial velocities at 0.296 s: (**a**) 0.2 m/s; (**b**) 0.3 m/s; (**c**) 0.4 m/s; and (**d**) 0.5 m

#### *4.2. The E*ff*ects of the Water Jet Flow Position*

Now, the liquid sloshing in a rectangular tank with a water jet flow from various positions at the top of the tank is compared. These correspond to the configurations shown in Figure 17, and there is also a harmonic rolling motion imposed on the geometry center of the tank with an amplitude 8◦ and a series of excitation frequencies, where P1 is the point of pressure monitoring. A1, A2, and A3 are three different locations for the water flow inlet, and d is 0.5, 0.3, and 0.1 m, respectively, which is the horizontal distance between the pressure probe and the water jet flow. The initial velocity of the water jet flow is 0.3 m/s.

**Figure 17.** Schematic view of the sloshing tank with water jet flows from different positions.

Figure 18 shows the comparisons of free surface profiles at *t* = 0.296 s with the water jet from three positions at the top of the tank when the excitation period is 1.85 s. From the contrast of the three free surfaces, it can be seen that the deformations of the free surface generated by the impact of three water flows are in a basic agreement, which generates the air cavities with the same shapes and sizes.

**Figure 18.** Comparison of free surface profiles at the instant *t* = 0.296 s with a water jet from three positions: (**a**) no water jet; (**b**) A1; (**c**) A2; (**d**) A3.

Figure 19 gives the peak value, which is obtained from the impact pressure of P1 at different excitation frequencies. It can be seen that the maximum pressure at P1 decreases when a water flow enters the sloshing tank with an initial velocity, and the excitation frequency where the maximum pressure occurs increases. Also, they change accordingly with the position of the water flow. When the horizontal distance between the water flow and P1 is closer, the effect is more obvious.

**Figure 19.** Comparisons of impact pressure with the water jet flow in different positions.

Through quantitative analysis, the maximum pressure is reduced by about 12.83% compared to the system with and without a water jet flow at A1, and this difference reaches up to 17.22% when the water jet flow comes from A3. It shows that the water jet flow can reduce the sloshing impact load on the tank wall. When the position of the water jet flow is closer to the pressure probe, the sloshing impact load becomes smaller. In addition, the water jet flow can also change the frequency where the maximum pressure occurs, so it is an effective way of avoiding the appearance of maximum pressure by using a water jet flow in an appropriate position.

#### *4.3. The E*ff*ects of the Water Jet Flow Number*

In this section, the effects of the water jet flow number are investigated based on the configuration shown in Figure 20. The same rectangular tank has a length of *L* = 1.2 m and a height of *h* = 0.6 m, which is partially filled up to an initial depth of 0.12 m. B1, B2, and B3 are three different locations for the water jet flow inlet. Again, a harmonic rolling motion is imposed on the tank with the same amplitude previously used. The sensor P1 is placed, as shown in the figure, to monitor the pressure variations. The same initial velocity of the water jet used in the previous section is also adopted here, which is 0.3 m/s.

**Figure 20.** Schematic view of the sloshing tank with a water jet: (**a**) only one water flow; (**b**) two water flows at B1 and B2; (**c**) two water flows at B1 and B3; (**d**) three water flows at B1, B2, and B3.

The particle snapshots of the pressure contour for the four cases at *t* = 0.296 s are shown in Figure 21, which illustrates the pattern of the free surface when the water jet flows enter the liquid. It shows that the water jet flow has played an important role in the pressure distribution of the sloshing process because it causes strong collisions between the water jet flow and the free surface below. It should be noted that there are some differences in the pressure distribution at different jet regions.

**Figure 21.** Comparison of free surface profiles at *t* = 0.296 s with a different number of water jets: (**a**) only one water flow; (**b**) two water flows at B1 and B2; (**c**) two water flows at B1 and B3; (**d**) three water flows at B1, B2, and B3.

Figure 22 shows the peak value of impact pressure at P1 obtained by ISPH at different excitation frequencies. The comparison includes the case of two water jet flows and only one water jet flow. In all the cases, the excitation frequency is a non-dimensional one, and its value is obtained between 0.9 and 1.5. All the cases correspond to the configurations shown in Figure 20a–c. By comparing the two water flows with only one water flow at B1, it can be seen that the maximum pressure at P1 can get a smaller value when another water flow is added at B2 or B3; meanwhile, the excitation frequency where the maximum pressure occurs also increases. Furthermore, the comparisons between different combinations of two water jet flows are shown in Figure 20b, c. The maximum pressure value of this case combined by B1 and B3 is the smallest. It is reduced by around 5.6% compared with only one water jet flow at B3 and is reduced by around 2.7% compared with the case combining B1 and B2. The excitation frequency where the maximum pressure occurs moves right again when B3 replaces B2. In summary, it demonstrates that the maximum pressure at P1 can be reduced when the other water jet flow is added, and it decreases more remarkably when the water flow is added closer to the pressure probe.

Figure 23 gives the peak value comparisons at P1 for the case of two water jet flows and three water jet flows. The excitation frequency adopts a non-dimensional one, which is between 0.9 and 1.5. The results are obtained according to the configurations shown in Figure 20b–d. The contrast demonstrates that it can obviously be improved for reducing the maximum pressure when more water jet flow is added.

**Figure 22.** Comparisons of impact pressure with the single water jet flow and two water jet flows.

**Figure 23.** Comparisons of impact pressure with two water flows and three water flows.

#### **5. Conclusions**

In this paper, the incompressible SPH method was used to simulate liquid sloshing in a tank with different configurations of the water jet flow, including the water jet flow with different initial velocities, positions, and water jet flow numbers. ISPH shows good agreement in both the free surface profiles and the impact pressures with the experimental data and demonstrates its great potential in predicting violent sloshing flows. The main purpose of this paper was to study the practical importance of the effect of a water jet flow on liquid sloshing through follow-on model applications of different configurations.

The main conclusions of the paper lie in the following aspects. Firstly, adding a water jet flow at the top of the sloshing tank can reduce the maximum impact pressure effectively. Secondly, adding a water jet flow at the top of the sloshing tank can change the excitation frequency where the maximum pressure occurs. Thirdly, adding different numbers of water jet flows also can decrease the value of maximum impact pressure. Finally, when the horizontal distance between water flow and pressure probe is closer, the effect on maximum pressure and excitation frequency is more obvious.

Furthermore, it should be noted that there are some limitations in the present sloshing model as follows:


**Author Contributions:** H.J. and Y.Y. made the computations and data analysis; Z.H. made the data analysis and did the proofreading; X.Z. did the proofreading and editing; Q.M. guided the engineering project and provided the data; Y.Y. drafted the manuscript with others. All authors contributed to the work.

**Funding:** This research work was fundedby the National Natural Science Foundation of China (Nos. 51879051; 51739001; 51579056, and 51639004); Natural Science Foundation of Heilongjiang Province in China (E2018024); Foundational Research Funds for the Central Universities (Nos. HEUCF170104; HEUCDZ1202); Defense Pre Research Funds Program (No. 9140A14020712CB01158).

**Acknowledgments:** The fifth author also thanks the Chang Jiang Visiting Chair Professorship scheme of the Chinese Ministry of Education, hosted by HEU.

**Conflicts of Interest:** The authors declare no conflict of interest.

#### **References**


© 2019 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (http://creativecommons.org/licenses/by/4.0/).

### *Article* **Simulating the Overtopping Failure of Homogeneous Embankment by a Double-Point Two-Phase MPM**

#### **Yong-Sen Yang, Ting-Ting Yang, Liu-Chao Qiu \* and Yu Han**

College of Water Resources & Civil Engineering, China Agricultural University, Beijing 100083, China **\*** Correspondence: qiuliuchao@cau.edu.cn; Tel.: +86-10-6273-6423

Received: 25 May 2019; Accepted: 5 August 2019; Published: 8 August 2019

**Abstract:** Embankments are usually constructed along rivers as a defense structure against flooding. Overtopping failure can cause devastating and fatal consequences to life and property of surrounding areas. This motivates researchers to study the formation, propagation, and destructive consequences of such hazards in risk analysis of hydraulic engineering. This paper reports a numerical simulation of failure processes in homogeneous embankments due to flow overtopping. The employed numerical approach is based on a double-point two-phase material point method (MPM) considering water–soil interaction and seepage effects. The simulated results are compared to available laboratory experiments in the literature. It was shown that the proposed method can predict the overtopping failure process of embankments with good accuracy. Furthermore, the effects of the cohesion, internal fiction angle, initial porosity, and maximum porosity of soil on the embankment failure are investigated.

**Keywords:** embankments; overtopping failure; material point method; water–soil interactions; numerical simulation

#### **1. Introduction**

Embankments are frequently used along rivers against flooding events, the probability of which has increased due to increasing intensity extreme weather events in recent years [1]. Their overtopping failure may cause serious loss of life and property downstream. For example, the Banqiao dam-breaking flood due to Typhoon Nina in 1975 lead to about 62,000 deaths [2]. Therefore, a deep investigation of overtopping failure mechanisms and processes in embankments is important for the risk management and assessment. Among the available study methods, numerical modeling has distinctive advantages such as the relatively low cost and the flexible application to almost arbitrary site configurations. In the past decades, many numerical investigations of the overtopping failure in embankments have been performed. For example, Powledge et al. [3] studied the embankment erosion due to overflow, Tingsanchali and Chinnarasri [4] used a one dimensional model to study the overtopping failure in embankment. Chinnarasri et al. [5] investigated the overtopping-induced progressive damage in embankment. Leopardi et al. [6] reviewed the modeling of embankment overtopping, Pontillo et al. [7] employed a two-phase model to evaluate the dike erosion induced by overtopping, Volz et al. [8] modeled the overtopping failure of a non-cohesive embankment by the dual-mesh approach, Mizutani et al. [9] simulated the overtopping failure of embankment considering infiltration effects, Guan et al. [10] and Kakinuma and Shimizu [11] developed a 2D shallow-water model for embankment breach, Evangelista [12] simulated dike erosion induced by dam-breaking flows using a depth-integrated two-phase model, Larese et al. [2] modeled the overtopping and failure of rockfill dams using the particle finite element method (PFEM). Although many efforts have been made to numerically reproduce the main features of the process, the overtopping failure of embankments is still poorly understood. Moreover, most of the aforementioned numerical methods are mesh-based approaches and have difficulties in accurately modeling the time evolution of the embankment breach, since the real overtopping failure phenomenon is, in fact, an unsteady process involving large

deformations, free-surface flows, moving boundary, and water–soil interactions. To deal with these challenges, several meshless methods, such as smoothed particle hydrodynamics (SPH) [13,14], moving particle semi-implicit method (MPS) [15], element-free Galerkin method (EFGM) [16], and material point method (MPM) [17,18], have been developed to model large deformation problems. For examples, Gotoh et al. [19,20] modeled embankment erosion due to overflow using MPS, Li et al. [21] simulated erosion of HPTRM levee based on SPH, Zhang et al. [22] modeled failures of dike due to water level-up and rainfall using SPH, Liu et al. [23] simulated piping erosion process of dike foundation by EFGM. Nikolic et al. [24] presented a discrete beam lattice model capable of simulating localized failure in a heterogeneous fluid-saturated poro-plastic solid. Zhao and Liang [25] applied MPM to model seepage flow through embankments, Martinelli et al. [26] modeled the failure of a sand dike due to seepage flow using two-phase MPM. In this work, the MPM was employed for modeling the overtopping failure of embankments because it has several advantages over other meshless methods such as EFGM, MPS, and SPH: The boundary condition realization is as easy as in the finite element method (FEM) for the use of background mesh in MPM [27]. Time-consuming neighbor particle searching is required in SPH, MPS and EFGM, but the MPM requires only the identification of particles relative to the background mesh. The MPM also avoids the tensile instability that is evident in the SPH [27].

The MPM was first proposed by Sulsky et al. [17] in 1994, which uses the Lagrangian particle and the Euler background mesh to describe large deformation problems. The particles move freely across the background mesh and carry all physical parameters such as stress, density, mass, velocity, and other historical variables. The MPM combines the advantages of the Euler method and the Lagrangian method to avoid the interference of convection terms and mesh distortion. It also can treat contact problems between different objects without any additional interface elements. In addition, the traditional soil constitutive laws such as Mohr–Coulomb, Drucker–Prager, and Cam Clay models can be easily implemented in MPM at the particle level [17,18]. Recent developments in the multiphase MPM formulations considering the soil–water interaction can be categorized into two main groups: Single-point approach and double-point approach. In the first approach, each material point (MP) possesses both the soil and water information. The water is described in the Eulerian framework while the soil skeleton is described in the Lagrangian framework. The single-point approach does not guarantee the mass conservation of water. In the second approach, two sets of MPs are used to represent water and solid maters. The double-point approach can guarantee both the solid and water mass conservation. Mackenzie-Helnwein et al. [28] used the double-point approach to model soil–water mixtures through a drag model for capturing the interphase interaction. Abe et al. [29] employed the double-point approach to model saturated soil without considering the relative acceleration between water and the soil skeleton. Bandara and Soga [30] developed a double-point two-phase MPM formulation based on the mixture theory considering the relative acceleration between water and the soil skeleton. This approach not only can be used for modeling rapid deformation problems but also guarantees both the soil skeleton and water mass conservation. It is therefore very suitable for modeling soil–water interaction problems such as currently considered embankment failure due to erosion.

The present work aims to numerically investigate the failure mechanism and failure process of homogeneous embankment due to overtopping flows based on a double-point two-phase MPM method.

#### **2. Mathematical Formulation**

The double-point two-phase MPM formulation [26,31] based on the mixture theory is used to model the water–soil interactions. In this approach, the soil skeleton and the pore liquid in the saturated soil are represented separately by two sets of material points: Solid material points and water material points (Figure 1). This paper focuses on the application of two-phase MPM to simulate embankment overtopping, only some of the basic equations solved are quoted for the convenience of readers in understanding the two-phase MPM formulation, and for detailed description of the algorithms readers are referred to the literature [26,31].

**Figure 1.** The schematic model for the double-point two-phase material point method (MPM) [32].

The governing equations of the double-point two-phase MPM include the mass and momentum conservations of both solid and liquid phases with their respective constitutive equations. The mass conservation equations for the solid and liquid phase are respectively given by [31].

$$\frac{D\mathbf{n\_L}}{Dt} = \mathbf{n\_S} \nabla \cdot \mathbf{v\_S} \tag{1}$$

$$\frac{D\varepsilon\_{\rm L}}{Dt} = \frac{1}{n\_{\rm L}} \left[ n\_{\rm S} \nabla \cdot \mathbf{v}\_{\rm S} + n\_{\rm L} \nabla \cdot \mathbf{v}\_{\rm L} + (\mathbf{v}\_{\rm L} - \mathbf{v}\_{\rm S}) \nabla n\_{\rm L} \right] \tag{2}$$

where *n*<sup>S</sup> is the volumetric concentration ratio of solid, *n*<sup>L</sup> is the volumetric concentration ratio of liquid. For saturated soils, *n*<sup>L</sup> is equivalent to porosity of the soil skeleton and *n*<sup>L</sup> + *n*<sup>S</sup> = 1. *n*<sup>L</sup> is a discontinuous function at the transition between free surface water and porous medium and it is characterized by two constant values: One in the free surface water and the value of the porosity *n* in the porous medium. Martinelli [33] proposed the concept of a transition zone in which an interpolated liquid concentration ratio is introduced which gives a smooth transition of *n*<sup>L</sup> between the free surface water and the porous medium. **v**<sup>S</sup> is the velocity vector of the solid phase, **v**<sup>L</sup> is the velocity vector of the liquid phase. ε<sup>L</sup> is the volumetric strain of the liquid. *<sup>D</sup>*(•) *Dt* denotes the material time derivative. Equation (1) represents the variation of volumetric concentration ratio of the liquid phase (i.e., porosity), and Equation (2) is also known as the storage equation and represents the volumetric strain rate of the pore liquid.

The momentum conservation of the solid and liquid phase can be respectively expressed as [31]:

$$
\overline{\rho}\_{\text{S}} \mathbf{a}\_{\text{S}} = \nabla \cdot \mathbf{\overline{\sigma}}\_{\text{S}} + \mathbf{f}\_{\text{L}}^{\text{d}} + \overline{\rho}\_{\text{S}} \mathbf{g} \tag{3}
$$

$$
\overline{\rho}\_{\rm L} \mathbf{a}\_{\rm L} = \nabla \cdot \overline{\mathbf{o}\_{\rm L}} - \mathbf{f}\_{\rm L}^{\rm d} + \overline{\rho}\_{\rm L} \mathbf{g} \tag{4}
$$

where ρ<sup>S</sup> and ρ<sup>L</sup> are respectively the partial densities of the solid and liquid, computed as the ratio of the mass of each constituent with respect to the reference volume. **g** is the gravity vector, **a**<sup>L</sup> is the liquid phase acceleration, **<sup>a</sup>**<sup>S</sup> is the solid phase acceleration. **¯** σ<sup>S</sup> = σ <sup>−</sup> *<sup>n</sup>*Sσ<sup>L</sup> and **¯** σ<sup>L</sup> = *n*Lσ<sup>L</sup> are the partial stresses of solid and liquid phases respectively, σ is the effective stress tensor, and σ<sup>L</sup> is the liquid phase stress tensor. **f** d <sup>L</sup> is the drag force, which represents the water–soil interaction due to the velocity difference between the two phases, it can be calculated by [34]

$$\mathbf{f}\_{\rm L}^{\rm d} = \frac{n\_{\rm L}^{2}\mu}{\kappa}(\mathbf{v}\_{\rm L} - \mathbf{v}\_{\rm S}) + \beta n\_{\rm L}^{3}\rho\_{\rm L}|\mathbf{v}\_{\rm L} - \mathbf{v}\_{\rm S}|(\mathbf{v}\_{\rm L} - \mathbf{v}\_{\rm S}) + \sigma\_{\rm L}\nabla n\_{\rm L} \tag{5}$$

where μ is the liquid dynamic viscosity, κ = *<sup>D</sup>*<sup>2</sup> *<sup>A</sup> <sup>n</sup>*<sup>3</sup> <sup>L</sup>/(1 − *n*L) <sup>2</sup> is the soil intrinsic permeability with *D* the average diameter of grains, β = *B*/ κ*An*<sup>3</sup> <sup>L</sup> is the non-Darcy flow coefficient, the empirical constants *A* and *B* are respectively set to 150 and 1.75 [34].

In order to fully describe the behavior of saturated soils, the constitutive equations for both phases are required. Assuming the validity of Terzaghi's effective stress concept, the general form of stress–strain relationship for soil skeleton is given by

$$\frac{D\sigma'}{Dt} = \mathbf{D}\frac{D\varepsilon}{Dt} \tag{6}$$

where **D** is the tangent stiffness matrix, σ is the effective stress vector, and ε is the total strain. For the liquid phase, the volumetric stresses is updated by

$$\frac{D\sigma\_{\text{Vol, L}}}{Dt} = K\_{\text{L}} \frac{D\varepsilon\_{\text{Vol, L}}}{Dt} \tag{7}$$

where *K*<sup>L</sup> is the liquid bulk modulus, εVol, L is the volumetric strain of the liquid. In addition, for pure liquid or fluidized mixture, the deviatoric part of the stress tensor is calculated by

$$
\sigma\_{\rm dev,L} = 2\mu \frac{D\varepsilon\_{\rm Vol,L}}{Dt} \tag{8}
$$

In the proposed double-point two-phase MPM, the saturated soils can be considered as a solid-like or liquid-like state according to the porosity. As shown in Figure 2a, in the case of low porosity, the solid skeleton grains are in contact and the soil behaves like a solid, its response thus can be modeled by a soil constitutive model. For a high-porosity soil, as shown in Figure 2b, the grains are not in contact and float together with the liquid phase and the soil response is modeled by the Navier–Stokes equation. In the present work, a pre-given maximum porosity *n*max is used to distinguish the two aforementioned states. When the soil porosity is less than the maximum porosity (*n*<sup>L</sup> = 1 − *n*<sup>S</sup> < *n*max), the mean effective stress decreases as the porosity increases and the mean effective stress vanishes once the grains are not in contact. When the porosity is larger than *n*max, fluidization occurs.

(**a**) solid-like response (low porosity) (**b**) liquid-like response (high porosity) **Figure 2.** Schematic diagram for soil behavior [34].

#### **3. Numerical Examples**

In this section, the proposed double-point two-phase MPM is first validated by an example of flow through a porous block. Then a numerical investigation of failure process in a homogeneous embankment due to flow overtopping was performed and the effects of the cohesion, internal fiction angle, initial porosity, and maximum porosity of soils on the embankment failure are investigated. All calculations are conducted with the software Anura3D\_v2016 (www.anura3d.com), which implements the double-point two-phase MPM formulation.

#### *3.1. Flow through a Porous Block*

In order to validate the capability of the current MPM algorithm to model the water–soil interaction, a numerical simulation of the flow through a porous block was performed and the calculated results were compared with the experimental data obtained by Liu et al. [35]. The experiment as sketched in Figure 3 was carried out in a transparent water tank with a size of 0.892 m × 0.44 m × 0.58 m, the porous block, which was made of gravel, was located at the center of the water tank (*x* = 0.30–0.59 m), and it was 0.29 m long, 0.44 m wide, and 0.37 m high. The water on the left side of the water tank was separated from the porous block by a movable gate with a thickness of 0.02 m, and the height was set to 0.35 m, and the porous block was confined in the initial region to ensure that the porous media was not allowed to move. The material parameters used in this simulation are listed in Table 1. The computational domain was divided into 4104 tetrahedral elements. Each mesh element in the porous block part contained four material points, while each element in the water contained 10 material points.

**Figure 3.** Geometry and discretization of the problem of flow through the porous block.


**Table 1.** The material parameters for flow through the porous block.

The simulation process was divided into two stages, the first stage was the initialization of stress, also called the gravity loading process. At this stage, the lateral displacement of the water was constrained and the stress initialization was done by increasing the gravity until the solution converged to a quasi-static equilibrium, which was checked by evaluating the normalized kinetic energy and the normalized out-of-balance force falls below a predefined tolerance as given in [31]. In addition, a local damping coefficient of 0.75 was applied to accelerate the convergence to quasi-static equilibrium

allowing a considerable reduction in the computational time. In the second stage, the lateral constraint of water was removed to simulate the instantaneous moving process of the gate, and then water flow impacted the porous block. Figure 4 shows some snapshots of the flow through the porous block at time of 0.35, 0.75, 1.15, and 1.95 s and the calculated free-surface profiles agree well with that of the experimental results [35], and comparisons of free surface profiles are shown in Figure 5. It is proven that the double-point MPM can simulate the process of flow through the porous medium considering the water–soil interaction.

**Figure 4.** Some snapshots of the flow through the porous block.

**Figure 5.** Comparison of free surface profile for flow through porous block: Simulation (-) and experiment (◦).

#### *3.2. Overtopping Failure of Homogeneous Embankments*

For a better understanding of the embankments failure mechanism and process due to the overtopping flow, in this section, MPM modeling of the overtopping erosion of a homogeneous embankment considering free water surface flow, water–soil interaction, and seepage effects was performed and the simulated results were compared with the experimental data obtained by Evangelista et al. [12]. The effects of the cohesion, internal fiction angle, initial porosity, and maximum porosity of soil on the embankment failure were investigated. As sketched in Figure 6a, this experiment was carried out in a Perspex horizontal channel with a 0.40 m wide rectangular cross-section and transparent walls. A 3.00 m long tank with a moving gate was used to store water, in which *H* represents the water level. The gate was suddenly lifted to produce the dam-break flow, which propagates over the downstream channel bottom. The first segment of the rigid bed was 1.00 m long and had a mild slope (5%), the subsequent one was a 2.00 m long, 0.05 m thick horizontal segment, on which the embankment was built.

(**b**) Geometry and MPM discretization

**Figure 6.** Experimental setup [12] and MPM model for the embankment overtopping problem.

According to Evangelista et al. [12], the friction at the lateral interface of the channel has no significant wall effects and the embankment undergoes an almost plane process because the flow overtops the entire embankment width. For these reasons, in our simulation, the overtopping failure process was simplified approximately as a two-dimensional problem. In addition, the gate opening was considered instantaneous following the work of Lauber and Hager [36]. The computational domain was divided into a set of tetrahedral elements with the size of 0.05 m. Each mesh element in the embankment part contained four material points, while each element in the water contained 10 material points (see Figure 6b). The water was modeled as a weak compressible fluid, and the Mohr–Coulomb model was used as the constitutive model of the soil. The material parameters used in this simulation are listed in Table 2. Three different cases with initial water level of *H* = 0.2, 0.25 and 0.395 m were conducted.

For the first case (*H* = 0.2 m), Figure 7 shows several snapshots of the simulated water surface and the embankment profiles along the channel at different times of 0.35, 0.8, 1.4, 2.0, and 3.0 s. Comparisons between the numerically simulated results and the experimental ones [12] are illustrated in Figure 8. Specifically, the embankment (black color) profiles and the water surface (blue color) at times 1.4 and 3.0 s, respectively, after the gate removal are plotted: Solid lines represent the numerical results and dotted lines indicate the experimental ones. In this case, the initial water level in the water tank was

almost equal to the top elevation of embankment, thus the gate opening–induced dam break flow did not have enough energy to overtop the dike embankment, and therefore, it induced erosion only on the upstream slope of the embankment.


**Table 2.** The material parameters of the two phase materials.

(**e**) *t* = 3.0 s

**Figure 7.** Snapshots of the simulated water surface and the embankment profiles along the channel (*H* = 0.2 m).

For the second case (*H* = 0.25 m), Figure 9 shows several snapshots of the simulated water surface and the embankment profiles along the channel at different times of 0.35, 0.8, 1.4, 2.0, and 3.0 s. Comparisons between the numerically simulated results and the experimental ones [12] are illustrated in Figure 10. Specifically, the embankment (black color) profiles and the water surface (blue color) at times 1.4 and 3.0 s, respectively, after the gate removal are plotted: Solid lines represent the numerical results and dotted lines indicate the experimental ones. In this case, the initial water level in tank was higher than the top elevation of embankment, thus the gate opening–induced dam breaking flow overtopped the embankment and eroded it along its entire profile. The upstream slope erosion was

more prominent due to the higher potential energy of the dam breaking flow. This dam breaking flow also reached the crest of the embankment this time, progressively eroding it as it passed.

**Figure 8.** Comparisons of the free water surface and embankment profiles between the simulation and experiment (*H* = 0.2 m).

(**e**) *t* = 3.0 s

**Figure 9.** Snapshots of the simulated water surface and the embankment profiles along the channel (*H* = 0.25 m).

For the third case (*H* = 0.395 m), Figure 11 shows several snapshots of the simulated water surface and the embankment profiles along the channel at different times of 0.35, 0.8, 1.4, 2.0, and 3.0 s. Comparisons between the numerically simulated results and the experimental ones [12] are illustrated in Figure 12. Specifically, the embankment (black color) profiles and the water surface (blue color) at times 1.4 and 3.0 s, respectively, after the gate removal are plotted: Solid lines represent the numerical results and dotted lines indicate the experimental ones. In this case, the potential energy of the gate opening induced dam break flow was about two times that in the previous case, thus leading to more prominent erosion in the embankment. The erosion was quite homogeneous along the entire profile, and the embankment profile shape was rounded during the process. During the process, steep portions of the embankment slopes underwent a combination of surface erosion and sliding failures due to the loss of supporting material.

**Figure 10.** Comparisons of the free water surface and embankment profiles between the simulation and experiment (*H* = 0.25 m).

**Figure 11.** Snapshots of the simulated water surface and the embankment profiles along the channel (*H* = 0.395 m).

**Figure 12.** Comparisons of the free water surface and embankment profiles between the simulation and experiment (*H* = 0.395 m).

In general, the erosion process of the embankment was greatly affected by the upstream water level elevation. The higher the initial water level, the more serious was the rate of damage of the embankment. By analyzing the embankment erosion process at low, medium, and high water levels, it can be seen that when the water level is low or medium (*H* = 0.2, 0.25 m), the water flow coming to the embankment takes a long time, and the erosion is mainly concentrated on the upstream slope surface. When the water level is high (*H* = 0.395 m), the water flow has a larger energy and rises rapidly along the embankment and overtops the crest to the downstream. Under the erosion and seepage of the water flow, the entire embankment is damaged seriously. Overall, the numerical simulation can relatively reproduce the overtopping erosion phenomenon, which validates the capacity of the proposed double-point MPM to model the overtopping failure of the homogeneous embankments with considering free water surface flow, water–soil interaction, and seepage effects.

In the previous simulation, the homogeneous embankments were considered as non-cohesive. However, we must remember that many embankments were built from cohesive materials, which behave quite differently than non-cohesive materials. In order to study the effect of cohesion on the embankment erosion process, the previous third case (*H* = 0.395 m) was selected, the simulations were conducted with three cohesions of *c* = 0, 5, and 10 kPa, the other conditions remained unchanged. Simulated results are shown in Figure 13, the erosion rate of the upstream and downstream slopes was more serious when *c* = 0 kPa. However, when *c* = 5 and 10 kPa, the rate of embankment erosion was very small and almost negligible. It can be seen that the increase of soil cohesion can effectively improve the resistance. As the cohesion increases, the erosion caused by the water flow is weakened, but when increased to a certain extent, the effect is weakened.

The effect of the internal fiction angle on the embankment erosion process was also investigated. Three different internal fiction angles (ϕ = 30◦, 36◦, 41◦) were used in the case of the initial water level *H* = 0.395 m, and the other parameters were the same as in the previous simulations. The calculated embankment profiles for different internal fiction angles at *t* = 1.4 s and *t* = 3.0 s are shown in Figure 14. Results show that the larger the internal friction angle is, the weaker the erosion will be. It is the internal friction angle that influences the contact between soil particles. For the small internal friction angle, it has weak resistance to overtopping flow, whereas a large friction angle can retard the flow.

**Figure 13.** The effect of the cohesion on the embankment profiles (*H* = 0.395 m).

**Figure 14.** The effect of the internal fiction angle on the embankment profiles.

In order to investigate the effect of the initial porosity on the embankment erosion process, three different initial porosities (*n* = 0.2, 0.3, 0.4) were selected to model in the case of the initial water level *H* = 0.395 m, and the other parameters were also the same as in the previous simulations. The calculated embankment profiles for different porosities at *t* = 1.4 s and *t* = 3.0 s are shown in Figure 15. It can be seen that in the early stage (*t* = 1.4 s) of the embankment failure, the erosion was not very obvious, only surface erosion caused by the water flow occurred and there was no infiltration. However, the profiles of the embankments at *t* = 3.0 s were different. The larger the porosity of the soil, the more obvious the embankment erosion. It indicates that the water flows began to infiltrate into the embankment in the middle stage, and the embankment erosion was induced by the water flow and seepage, and the seepage failure mainly occurred in the middle and later stages. With the gradual infiltration of the water flow, the water flow infiltrated completely in the later stage, and the pores of the soil were filled with water, the soil became saturated. When the effective stress of the soil skeletons is 0, the soil may be liquefied and moved with the water flow. Hence, decreasing the porosity of the soil can slow down the infiltration rate and reduce the damage of the embankments.

**Figure 15.** The effect of the initial porosity on the embankment profiles.

In the proposed double-point two-phase MPM, a pre-given maximum porosity *n*max was used to distinguish the saturated soils as a solid-like or liquid-like state. When the porosity was larger than *n*max, fluidization occurred. The grains are not in contact and float together with the liquid phase. Here, three cases of different porosities (*n*max = 0.45, 0.50, 0.55) were conducted to study the effect of maximum porosity on the embankment erosion process, and the other conditions were the same as in the previous simulations. The calculated embankment profiles for different cases at *t* = 1.4 s and *t* = 3.0 s are shown in Figure 16. Results show that the maximum porosity plays an important role on the dam failure due to overtopping flow. The larger the maximum porosity is, the more difficult it is for soil particles to reach the fluidization state. Hence, it is very important to get the right maximum porosity in the simulations, and more accurate simulation results can only be obtained with the appropriate value.

**Figure 16.** The effect of maximum porosity on the embankment profiles.

#### **4. Conclusions**

In this paper, a double-point two-phase material point method was employed to model problems involving free water surface flow, large deformations, and water–soil interaction. The proposed method was first validated by the example of the flow through a porous block problem. Then a numerical investigation of failure process in homogeneous embankments due to overtopping flow was performed. By comparing the water flow and embankments profiles with the experimental data, it was shown that the current double-point two-phase MPM can predict the overtopping flow pattern and the development of embankment erosion with good accuracy, and the higher the initial water level, the more obviously and faster the embankment is eroded. Furthermore, the effects of the cohesion, internal fiction angle, initial porosity, and maximum porosity of soil on the embankment failure were investigated. The results show that these parameters have important influence on the erosion process. With the increase of the cohesion, the resistance of the embankment can be improved to a certain extent, but the effect of the cohesion is gradually weakened. The larger the internal friction angle is, the weaker the erosion will be. In addition, the larger the initial porosity of the soil, the more obvious the embankment erosion. The larger the maximum porosity is, the more difficult it is for soil particles to reach the fluidization state. These investigated results indicate that the double-point two-phase MPM is capable of predicting and reproducing the failure process of embankments due to the overtopping flow. The proposed method is an alternative promising tool for investigating complex failure mechanism in water–soil interaction.

**Author Contributions:** Y.-S.Y. and T.-T.Y. conducted the numerical computations and data analysis; Y.-S.Y. and L.-C.Q. drafted the manuscript; Y.H. did the proof reading and editing. All authors contributed to the work.

**Funding:** This research was funded by the National Key R & D Program of China (No. 2018YFC0406604) and the National Natural Science Foundation of China (No. 11772351).

**Acknowledgments:** The authors would like to acknowledge Anura3D MPM Research Community for providing the external version of the software Anura3D\_v2016.

**Conflicts of Interest:** The authors declare no conflict of interest.

#### **References**


© 2019 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (http://creativecommons.org/licenses/by/4.0/).

### *Article* **Local Scour of Armor Layer Processes around the Circular Pier in Non-Uniform Gravel Bed**

**Manish Pandey 1, Su-Chin Chen 2,\*, P. K. Sharma 3, C. S. P. Ojha <sup>3</sup> and V. Kumar <sup>4</sup>**


Received: 27 June 2019; Accepted: 8 July 2019; Published: 10 July 2019

**Abstract:** Flume experiments have been carried out under clear water scour conditions to analyze the maximum equilibrium scour depth and scour processes in armored streambeds. A total of 85 experiments have been carried out using different diameters of circular piers and non-uniform gravels. A graphical approach for dimensionless scour depth in equilibrium condition around the circular pier in armored streambeds has been developed. As per this curve, the maximum dimensionless scour depth variation with dimensionless armor particle size depends on the densimetric particle Froude number (Fr*<sup>d</sup>*<sup>50</sup> ), and the decreasing rate of dimensionless scour depth decreases with the value of Fr*<sup>d</sup>*<sup>50</sup> .

**Keywords:** pier scour; non-uniform sediment; armor layer; equilibrium scour depth processes; clear water scour condition

#### **1. Introduction**

Scour around bridge elements like abutments and piers have been recognized as the main cause of bridge collapse [1]. Hence, for the economical and safe design of bridge elements, it is necessary to estimate precise calculations of scour depths. Several studies are available on this subject, mostly dealing with uniform non-cohesive sediment beds [2–9]. For refraining bridges with pier hazards, it is necessary to propose a pier scour study based on many factors, classified as structural parameters of pier, bed sediment characteristics and hydraulic conditions prevailing at the site. So far, those objectives have not been met.

Characteristically, riverbeds in upper reaches deal with non-uniform coarser sediments such as non-uniform gravels, sand-gravel mixtures, sand-boulder-gravel mixtures etc. However, present study only focuses on non-uniform gravels. It is hard to accurately predict scour depths using uniform sediment bed approaches [10]. Sediment transport difficulties involving non-uniform gravel beds could not be simply explained critically. This is because the physical process of sediment transport phenomenon has not yet been ideally defined. Jueyi et al. [11] stated that a protecting layer of coarser particles formed around the pier is identified as the armor layer. In other words, the formation of an armor layer in a streambed stops the scouring process; however, some finer particles from the stream can be transported downstream through the armored bed [12]. An armor layer plays a significant role in river evolution.

When the incipient condition of uniform sediment in a streambed exceeds, the particles start to move. However, the entrainment development of a single particle in a non-uniform sediment bed is complex. The resistance to single particle motion is a function of shape, size and relative density of the particle. Schmidt and Gintz [13] stated that Platy shaped particles have a lesser possibility to be entrained than more compacted sediments. In addition, resistance to particle movement is also affected by the exposure and hiding effects in non-uniform sediments [13].

The development of scour processes in non-uniform streambeds stops because of the development of a steady armor coating, due to the removal of fine sediment particles by the flowing water from the channel bed material (parent bed material) near the pier. Usually, the value of geometric standard deviation of the stratified layer particles is less than parent bed material and these particles are larger in median diameter than parent bed particles [14]. It was observed experimentally by Jueyi et al. [11] that the removal of finer sediment tends to occur earlier because it entails less energy to move. The development of armor layer by flushing fine particles from the scoured hole, and gradual exposure of coarser particles does not allow the further removal of sediment from the scoured hole. This phenomenon is also known as the end condition for scour in non-uniform sediments, and the scoured zone does not undergo further changes with time. The condition of end scour and the formation of the armor layer is dependent on the approach flow parameters, sediment properties, channel bed geometry and pier shape.

The armor layer at equilibrium scour condition on streambed plays a significant role in pier scour. In the literature, only a few studies are available related to scour processes and the behavior of different parameters around bridge piers in armored streambeds [11,15–17]. Most recently, Jueyi et al. [11] derived a maximum dimensionless scour depth equation in the presence of the armoring layer. Guo [15] and Kim et al. [16] conducted pier scour study in non-uniform sediment and proposed maximum scour depth relationships. However, they did not consider the armor layer in their study. Jueyi et al. [11] conducted an experimental study for abutment scour in non-uniform gravel bed and proposed an empirical relationship to describe the maximum depth of scour at equilibrium state.

Presently, only a few studies are available to describe the scour processes around bridge piers in non-uniform streambeds [10,14–16]. In this study, we investigated the scour processes and behavior of different parameters on maximum scour depth in armored streambed. All experiments have been carried out in a laboratory flume under the clear water scour condition using different sizes of non-uniform gravels and flow conditions.

#### **2. Experimental Setup and Procedure**

In this study, a total of 85 tests were carried out in a 20.0 m long and 1.0 m wide rectangular flume, as shown in Figure 1. Four different diameters (*b*) of circular bridge piers, i.e., 6.6 cm, 8.4 cm, 11.5 cm and 13.5 cm were used for present experiments. All piers were made of hollow cast iron pipes.

**Figure 1.** Schematic view of the experimental setup.

All experiments have been completed under different flow conditions and sediment properties like time-averaged velocities (*U*), approach flow depths (*y*), pier diameters (*b*), median diameter of parent streambed particles (*d*50) and geometric standard deviation (σ). Test limits are mentioned in Table 1 for cases (a–h).


**Table 1.** Experimental data of non-uniform gravel bed.

The working section starts from 6.0 m of the flume entrance, as shown in Figure 1. The working section (8.0 m long and 0.35 m deep) was completely filled with eight different sizes of non-uniform gravels (σ > 1.4) having 3.26, 3.82, 4.38, 4.94, 7.5, 8.6, 9.1 and 10.7 mm median diameter (*d*50) and 1.57, 1.66, 1.88, 1.56, 1.59, 1.63, 1.71 and 1.51 geometric standard deviation of parent particle size distribution, respectively. The discharge was measured with an ultrasonic flow meter, which was provided at inlet pier of the flume. Approach flow depth was adjusted using a tailgate, which was located at the downstream end of the flume. A wave regulator or straighter was facilitated at the flume entrance to produce a uniform or near uniform flow condition in the experimental flume. In the present study, the maximum scour depth (*hs*) was measured with a Vernier point gauge. All tests were prepared for 24 hours. However, experimentally it was observed that the equilibrium scour stage was reached within 14 hours. After the equilibrium time of scour" i.e., 14 hours, the scour depth at different points near the pier was the same at every 30 minutes' interval. Before the start of each experiment, the working section of the flume was made perfectly level with respect to flume bed and covered with a thin Perspex sheet. Once preset flow conditions were achieved, the Perspex sheet was separated very carefully to avoid undesirable scour around the pier. Figure 2 shows a definition sketch of pier scour in armored bed for case (b), followed by Table 1.

**Figure 2.** Scour around circular pier in armored streambed for case (b), refer to Table 1.

After equilibrium scour state, the maximum scour depth (*hs*) was measured with a Vernier point gauge at upstream nose of the pier. We also collected the armor layer sediment for further particle size distribution and used the median diameter of armor layer (*d*50*a*) for scour depth analysis.

Critical velocity *Uc*, was computed by Melville and Sutherland [18] and varied with present flume experiments. Shield's approach was obtained for calculating the critical shear velocity (*u\*c*).

#### **3. Scour Processes and Armor Layer Formation**

The interaction between gravel particles and the flow is complex, and this a significant information for designing a pier in gravel-bed streams. Armoring is a common phenomenon in rivers and occurs when selective entrainment of finer particles leave the stream-bed with a coarser gravel size than that of the underlying bed. Experimentally, it was stated by various researchers that when shear velocity (*u\**) exceeds the critical shear velocity (*u\**c) of streambed particles then incipient motion condition develops in the stream, and then further increment in shear velocity develops scour processes near the bridge elements [6,14]. Melville and Sutherland [18] stated that the scour near the bridge elements occurs when approach shear velocity is almost 0.5 times of *u\*c* of armor bed particles. At the critical condition, the selective entrainment of gravels leaves the higher, more stable gravels on the stream-bed. Further increase in flow, the entrainment frequency increases and an armor layer of coarser gravel particles starts to form. The scour hole starts forming due to the formation of primary-vortex and drops into the scour hole as the scour hole volume increases with time [16]. Hydrodynamic forces are induced on particles and eliminate the particles around the pier. The drag and lift forces decrease with a rise in scour depth. In actuality, the depletion of the particles near the bridge elements (pier and abutments) occur layer-by-layer [14]. The scour processes end in non-uniform sediment beds with the development of a stable armor layer [16]. At equilibrium scour condition, this armor layer does not permit the further removal sediment from the scoured hole. An armor layer holds partially of cluster creations, which are more constant than specific particles. The cluster creations have a noteworthy influence on pier scour, they can offer more stable bed after equilibrium scour state. After flume experimental study, it has been found that the stable armor layer usually has uniform sediment, with lesser geometric standard deviation in comparison to parent streambed material and depends on the strength and quantity of coarser particles [14].

The strength of the scour hole depends on movement of coarser particles [11]. In this study, it has been noted by us that coarser particles are more stable inside the scour hole because of their larger mass (more energy is required to remove from scour hole). At the same phase, finer gravels are protected by coarser gravels. Slowly, this process results in the development of clusters of different sizes of gravels. For a definite expanse, these clusters show the same characteristics [11]. Kothyari et al. [14] stated that the formation of armor layer may be noticed as a nonstop process of developments and collapses of clusters.

#### **4. Results and Discussion**

The stability of armor layer around the pier is influenced by parent particle size, median diameter of armor particles and approach flow properties. The maximum scour depth (*hs*) variations at equilibrium condition with time-averaged flow velocity (*U*) for different diameters of piers along with different sizes of gravels is shown in Figure 3a–h. It can be noted that coarser gravels show less variation of maximum scour depth with approach flow velocity, while medium gravels show higher variation. It can also be noted from Figure 3, the maximum scour depth in armored streambeds increases with approach flow velocity and that for any sizes of streambed particles, different sizes of particles of armored streambeds illustrate a different maximum scour depth. It has been found that the maximum scour depth variation is higher for fine gravel beds (*d*<sup>50</sup> = 3.26 − 4.94 mm) and lower for coarser gravel bed (*d*<sup>50</sup> = 7.5 − 10.7 mm). For particular streambed sediment, the maximum scour depth increases with increase in pier diameter and approach velocity.

**Figure 3.** Maximum scour depth variation with approach flow velocity.

We observed that the median diameter of armoring layer (*d*50*a*) at equilibrium scour condition depends on σ and increases with σ, as shown in Figure 4. Armoring ratio (*d*50a/*d*50) defines the degree of armoring and mainly depends on the flow properties [19]. For σ = 1.63 − 1.88, the armoring ratio varied from 1.7–2.5, while for the lesser values of σ, i.e., σ = 1.5 − 1.6, the armoring ratio varied from 1.25–1.7, as can be seen in Figure 4. Wu et al. [17] stated that the size of equilibrium scour-hole illustrates a strong dependency with parent particle size distribution. Kothyari et al. [14] concluded that the maximum scour depth variation in non-uniform and non-cohesive sediment mixture is a function of geometric standard deviation of parent particle size distribution. Similar results were noted by us in the present study and the value of armoring ratio increases with σ. However, the increment in dimensionless scour depth with armor ratio is very gradual up to *hs*/*RL* = 0.5, where *RL* is the reference length and defined as (*b2y*) <sup>1</sup>/<sup>3</sup> by Oliveto and Hager [20,21] and Kothyari et al. [14]. After *hs*/*RL* = 0.5, there is no change in dimensionless scour depth with armor ratio, as can be seen in Figure 4. It indicates that the dimensionless scour depths have -less dependency on degree of armoring ratio, but maximum scour depths have shown noteworthy changes with pier diameter and time-averaged velocity, as can be seen in Figure 3.

**Figure 4.** Dimensionless scour depth variation with armor ratio.

Figure 5 shows the dimensionless scour depth (*hs*/*RL*) variation with *y*/*d*<sup>50</sup> ratio, where *y* is approach flow depth. For a particular sediment property, the maximum scour depth increases with approach flow depth (*y*). Figure 5 clearly shows that the maximum scour depth variation with approach flow depth is higher for *d*<sup>50</sup> = 3.26 − 4.94 mm (fine gravel) and lowest for *d*<sup>50</sup> = 7.5 − 10.7 mm (coarse gravel). Moreover, the finer streambed particles (*d*<sup>50</sup> = 3.26 − 4.94 mm) result in a higher dimensionless scour depth variation than that of coarser streambed particles (*d*<sup>50</sup> = 7.5 − 10.7 mm).

Figure 6 shows the dimensionless scour depth (*hs*/*RL*) variation with *d*50*a*/*RL*, for different ranges of critical velocity ratio (*U*/*Uc*), and critical velocity ratio is known as the ratio of the approach mean velocity to the critical velocity of sediment. The range of critical velocity ratio is from 0.70 to 0.98. It can be identified that the dimensionless scour depth variation with *d*50*a*/*RL* is highest for maximum range of critical velocity ratio, i.e., *U*/*Uc* ≈ 0.95 − 0.98. Hence, maximum number of experiments were carried out for *U*/*Uc* ≈ 0.90 − 0.95 and it was noted that for all ranges of critical velocities, *hs* increases with decrease of *d*50*a*. It has been observed that the lowest variation of dimensionless scour depth with *d*50*a*/*RL* for *U*/*Uc* ≈ 0.70 − 0.80, as can be seen in Figure 6. Figure 6 clearly illustrates that the dimensionless scour depth variation increases with *U*/*Uc*. It can also be noted that the time to reach near equilibrium depths for low values of *U*/*Uc* is larger than the higher values of *U*/*Uc*.

**Figure 5.** Dimensionless scour depth variation with *y*/*d*50.

**Figure 6.** Dimensionless scour depth (*hs*/*RL*) variation with *d*50*a*/*RL.*

#### *Implications for Design*

In the past studies, several investigators have derived numerous dimensionless scour depth relationships [2,14,22–28]. Present experimental results show that scour depth is not only a function of flow parameters, but also depends on particle size and geometric standard deviation of parent bed material and armor layer particle size. Pier scour experimental data collected under armor streambeds are shown in Table 1, and were used to illustrate the dimensionless scour depth variation with parent densimetric particle Froude number, Frd50 = *U*/[(*S* − 1)*gd*50] 0.5 and ratio of armor layer particle size − reference length (*d*50*a*/*RL*), as shown in Figure 7. Where *S* is specific gravity of gravel equal to 2.65. It can be noted that for a particular range of Frd50 , the dimensionless scour depth decreases with an increase in *d*50*a*/*RL*. The reducing frequency of dimensionless scour depth was very rapid for the higher value of densimetric Froude number, as can be seen in Figure 7. It can be noted that the decreasing rate of dimensionless scour depth variation with *d*50*a*/*RL* is lowest for Frd50 = 1.6–1.7. Meanwhile, Figure 5 shows the dimensionless scour depth gradually increasing with increase in armor ratio, for a particular sediment property. Figures 6 and 7clearly illustrate that the dimensionless scour depth variation with *d*50*a*/*RL* depends on *U*/*Uc* and Frd50 , respectively. During the experiments, it was observed that the armor layer around the pier is washed out for higher value of *U*/*Uc* due to the higher critical velocity of armor particles than the parent particles. For an armor-layered streambed, maximum scour depth increases with the decrease of parent streambed particle size and the increase of pier diameter and approach velocity, as shown in Figures 3 and 5. Therefore, under fixed flow conditions and parent streambed properties, the maximum scour depth in equilibrium condition is larger in finer parent particle streambeds.

**Figure 7.** Dimensionless scour depth variation with *d*50*a*/*RL* and Fr*d*50.

#### **5. Conclusions**

To date, very few studies are available for local scour around pier founded in armored gravel-bed streams. In the present study, a total of 85 experiments were carried out under different armor streambeds, along with different pier diameters and flow conditions. Using these different parameters, the maximum scour depth in armor streambeds was analyzed. A new graphical approach was proposed for the maximum scour depth computation.

It was observed that the scour hole around the pier occurs through the formation of an armor layer and scour hole has been reached within equilibrium scour state, when the armor layer particles cover the scour hole and no further movement allows from parent bed material. We found that the maximum scour depth is not only a function of flow parameters, but also depends on particle size and geometric standard deviation of parent bed material and armor layer particle size. The dimensionless scour depth variation with *d*50*a*/*RL* depends on the parent densimetric particle Froude number, and the decreasing rate of dimensionless scour depth decreases with Frd50 . We found that the influence of armor ratio on dimensionless scour depth is less, up to *hs*/*RL* = 0.5, however, after *hs*/*RL* = 0.5, change in dimensionless scour depth was negligible. The maximum scour depth around the pier increases with increase in pier diameter, approach velocity, critical velocity ratio and decreases with the size of parent streambed particles. The maximum scour depth decreases with an increase of armor layer particle size for a particular bed material having the same parent bed material. This research indicates the necessity for additional non-uniform coarser bed scour study as it relates to fluvial hydraulics.

**Author Contributions:** M.P., P.K.S. and C.S.P.O. conducted a series of experiments; M.P., V.K. and S.-C.C. analyzed data, applied methodology and validated this experimental study; S.-C.C., P.K.S., V.K., M.P. and C.S.P.O. reviewed and edited; funding acquisition by S.-C.C.

**Funding:** This research was funded by the Ministry of Education under the Higher Education Sprout Project.

**Acknowledgments:** The authors kindly acknowledge the staff of hydraulic laboratory, I.I.T. Roorkee and students of department of soil and water conservation, National Chung Hsing University, Taichung, Taiwan.

**Conflicts of Interest:** The authors declare no conflict of interest.

#### **Abbreviations**


#### **References**


© 2019 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (http://creativecommons.org/licenses/by/4.0/).
