Abstract

Sediment scouring is a common example of highly dynamic sediment transport. Considering its complexities, the accurate prediction of such a highly dynamic multiphase granular flow system is a challenge for the traditional numerical techniques that rely on a mesh system. The mesh-free particle methods are a newer generation of numerical techniques with an inherent ability to deal with the deformations and fragmentations of a multiphase continuum. This study aims at developing and evaluating a multiphase mesh-free particle model based on the weakly compressible moving particle semi-implicit formulation for simulation of sediment scouring. The sediment material is considered as a non-Newtonian viscoplastic fluid, whose behavior is predicted using a regularized μ(I) rheological model in combination with pressure-dependent yield criteria. The model is first validated for a benchmark problem of viscoplastic Poiseuille flow. It is then applied and evaluated for the study of two classical sediment scouring cases. The results show that the high-velocity flow currents and the circulations can create a low-viscosity region on the surface of the sediment continuum. Comparing the numerical results with the experimental measurements shows a good accuracy in prediction of the sediment profile, especially the shape and dimensions of the scour hole.

INTRODUCTION

Highly dynamic (or rapid) movement of sediment material is a problem of common occurrence in many hydro-environmental and engineering settings. Reservoir flushing, dam break flow over sediment bed, sediment scouring around structures (e.g., under pipes or around bridge piers), riverbank failure, and underwater landslides, are a few examples of the highly dynamic sediment movement, characterized by rapid motion of dense sediment mass with large deformations and fragmentation and frequent topology change. Among these examples, scouring has a special importance due to its effect on the stability of the hydraulic structures. The highly erosive flow (with high velocity and strong flow circulation) is often the driver of sediment scouring. The foundations of contemporary research on sediment processes (e.g., erosion, transport, deposition), can be found in many experimental, field, and soft-computing studies. Sediment scouring has also been well studied using these techniques (e.g., Azamathulla & Ghani 2010; Dogan & Arisoy 2015; Guan et al. 2016; Park et al. 2016; Roushangar et al. 2017; Barati et al. 2018).

With the advancement in computing powers and the numerical algorithms, the numerical simulations have gained much attention. Traditional approaches for simulation of sediment processes have mostly been based on single-phase modeling, in which the flow equations are coupled with the advection-diffusion equation (for suspended sediments) and the empirical or semi-empirical relationships (van Rijn 1984) for determining processes such as bed-load transport. These methods, which are the base of the majority of the existing models, have proved their effectiveness in some past studies (e.g., Falconer & Owens 1990; Kashyap et al. 2017).

Nevertheless, they are not suitable for dealing with the highly dynamic movement of sediments (the case of sediment scouring) as they are not able to predict the underlying physics in detail. Multiphase models, with separate water and sediment phases, have proven to be an enabler for many sediment transport simulations. These methods are based either on continuum-based mesh-based Eulerian or Lagrangian methods such as finite element method, finite volume method, and material point method (Hsu et al. 2003; Larese et al. 2013; Lee et al. 2016), or discrete-based Lagrangian methods such as discrete element method (Calantoni et al. 2004; Tao & Li 2018). The first group is ill suited in dealing with the multiphase flows with large interfacial deformations and fragmentations due to mesh-adaptability and connectivity problems (Shakibaeinia & Jin 2011) and the second group is computationally expensive for large-scale problems as they simulate each individual grain of sediment (Rycroft et al. 2010). The mesh-free particle (Lagrangian) methods developed for continuum mechanics problems, such as smoothed particle hydrodynamics (SPH) (Gingold & Monaghan 1977) and moving particle semi-implicit method (MPS) (Koshizuka & Oka 1996), provide simple and natural yet powerful tools for the simulation of highly dynamic multiphase flows. The mesh-free particle methods can accurately capture any interfacial deformations and fragmentations. Therefore, they can be ideal candidates for dealing with the complexities involved in highly dynamic sediment movements. They have the Lagrangian nature of discrete-based methods, yet do not have their scalability issue, as they model the multiphase continuum.

Particularly, the MPS, on which this study is based, has proven its capability for many complex fluid flow problems (e.g., Gotoh & Sakai 2006; Ataie-Ashtiani & Farhadi 2006; Shakibaeinia & Jin 2012a, 2012b; Gambaruto 2015). Shakibaeinia & Jin (2010) proposed a weakly compressible version of the MPS method (WC-MPS) to improve efficiency and also the numerical fluctuations associated with the original MPS. Jafari-Nodoushan et al. (2016) improved the accuracy of the model near the solid and inflow boundaries. Shakibaeinia & Jin (2012a) proposed a multiphase MPS method for multi-density multi-viscosity systems. Fadafan & Kermani (2017) improved pressures' stability properties with the MPS method. MPS is also very versatile in the simulation of non-Newtonian flows, by adapting various rheological relations. By combining the multiphase MPS model with visco-plastic rheological relations such as Herschel–Bulkley and Bingham plastic, Shakibaeinia & Jin (2011, 2012b) dveloped a model for simulation of flow-driven sediment flows for the example application of mobile-bed dam break and sand jets and plumes. Fu & Jin (2015) and Nabian & Farhadi (2016) used a similar approach for simulation gravity-driven sediment flow for the case of submerged landslides. Similarly, the SPH method was used for simulation gravity-driven (Capone et al. 2010) and flow-driven sediment flow (Manenti et al. 2011; Khanpour et al. 2016). More recently, Jafari-Nodoushan et al. (2018) and Tajnesaie et al. (2018) have developed and fully characterized a multiphase WC-MPS model for gravity-driven sediment flows. Unlike past studies, they were based on a pressure-dependent viscoplastic model with a dynamic inter-grain mechanical pressure.

Another important aspect of the continuum-based numerical models is the constitutive (rheological) model used to mathematically describe the non-Newtonian behavior of a granular (sediment) continuum. While the general visco-plastic models such as Bingham plastic, and Herschel–Bulkley have been two of the widely used models in past studies, more recent attention has been directed towards μ(I) rheological model (Jop et al. 2006). Unlike the Bingham plastic and Herschel–Bulkley models, the μ(I) rheological model is specialized for granular flows. In this model, the rheological properties are explicitly related to the the physical properties of granular material (e.g., grain size and density), and the granular flow regimes (quasi-static, inertial, and viscous regimes) are taken into account. Furthermore, Tajnesaie et al. (2018) showed that the μ(I) model produces slightly more accurate results than the Herschel–Bulkley model.

The present research therefore aims at developing and implementing a mesh-free particle (Lagrangian) model based on multiphase WC-MPS formulation for the simulation of two-phase highly dynamic movement of sediment for the specific case of sediment scouring. The method treats the multiphase system as a single multi-viscosity and multi-density continuum. An exponentially regularized μ(I) rheological model with the fully dynamic inter-grain mechanical pressure (unlike past studies) is applied to predict the non-Newtonian behavior of the sediment phase. The benchmark case of viscoplastic Poiseuille flow is used to evaluate the capabilities of the model dealing with non-Newtonian viscoplastic behaviors. The model is then validated and evaluated for simulation of the sediment scour around a submerged pipe and scour resulting from water jets (in a plunge pool) with different angles.

GOVERNING EQUATIONS

The governing equations for a multiphase flow system in a Lagrangian frame, including mass and momentum conservation and motion in Lagrangian coordinates, are as follows: 
formula
(1)
where, t is time, r = r(x, y) is the position vector, u = u(u, v) is velocity vector, ρ is density, g is gravity acceleration, p is pressure, and τ is deviatoric stress (shear stress) tensor given by: 
formula
(2)
where I is the unit tensor, η is the effective viscosity of materials (given by the rheological model), ξ is the second coefficient of viscosity, and E is the rate of strain tensor. Here, in this study, the flow is considered weakly compressible and an equation of state (EOS) predicts the pressure using the density field p=f(ρ). For a weakly compressible fluid, the divergence of the velocity vector goes to zero , therefore the stress tensor reduces to τ =2ηE. Considering the two-phase system of sediments and water as a single multi-density multi-viscosity continuum, one set of governing equations (Equation (1)) is used for both phases (similar to Shakibaeinia & Jin (2012a) and Hu & Adams (2007)) and no additional interfacial forces will be required.

Sediment rheology

To calculate the stress tensor, one needs to calculate the effective viscosity η. While the effective viscosity of water phase is constant (independent of time and stress) that of sediment phase is determined through a rheological model (constitutive relationship). Among the rheological models and theories developed for granular and sediment flows, a frictional model based on the μ(I) law (Jop et al. 2006; Forterre & Pouliquen 2008) has recently gained wide attention. The μ(I) rheological model is derived from the analysis of a large number of laboratory and numerical data, providing a general framework for describing a wide range of granular material behavior. Before defining the μ(I) model, one should remember the required basic definitions of tensor calculus, in which for D-dimensional arbitrary tensor A =[akl], the trace is given by , transposed by AT, the first invariant by AI, second invariant by AII, and magnitude by . The μ(I) model is a viscoplastic model, in which the material behaves as a rigid body for stresses less than a yield stress, τy and flows as a viscous fluid for stresses larger than τy. This model defines the effective viscosity as: 
formula
(3)
where μ(I) is the inertial-dependent coefficient of friction, and ps is mechanical pressure (normal stress) between the sediment grains. Here, the yield stress, τy is calculated from the Drucker–Prager yield criterion (Chen & Ling 1996), which for cohesionless material is calculated as τy = μ1ps where μ1 represents the first friction coefficient (a tangent of angle of repose θr). The coefficient of friction μ(I) depends on a dimensionless number, called the inertia number I, defined as the ratio between microscopic and macroscopic time scales, and depending on the regimes of granular flow. The coefficient of friction is given by (Jop et al. 2006): 
formula
(4)
where μ2 is the second friction coefficient (μ1 and μ2 represent the upper and lower limits of μ(I)), I0 is a reference inertial number (a constant coefficient, commonly determined through experiments). For a submerged granular flow with a low-viscosity ambient fluid (Stokes number = ds(ρsps)0.5/ηf >> 1) the inertial number is given by (Forterre & Pouliquen 2008): 
formula
(5)
where ds is sediment grain size and cd is the drag coefficient, and ρf and ηf are the density and viscosity of the ambient fluid, respectively. Therefore, the effective viscosity can be calculated by: 
formula
(6)
The equation is very similar to the generalized visco-plastic rheological relations such as Herschel–Bulkley and Bingham plastic (see Shakibaeinia & Jin 2011; Tajnesaie et al. 2018). Nevertheless, it relates the post-failure rheology to a variable term, which not only is a function of the material properties, but also depends on the pressure and strain rate. Equation (6) is singular and discontinues when the strain rate goes to zero (E → 0), therefore it cannot be directly used in the numerical model. To avoid this issue, a regularized continuous form has to be used. Here, we combine it with a widely used exponential regularization proposed by Papanastasiou (1987) to Equation (7) as: 
formula
(7)
where m is an exponent that controls the viscosity growth. This equation is valid for stresses smaller and greater than the yields stress and is continuous as for E → 0 the effective viscosity goes to the finite maximum value of .

NUMERICAL METHOD

The mesh-free numerical method of this study is based on the moving particle semi-implicit (MPS) formulation. In MPS, the computational domain is represented and discretized with a set of free-to-move computational nodes, called particles (Figure 1). The field variable is assigned to each particle. Interpolation of the quantities and integration of the spatial derivatives on each particle is through a weight averaging process between the target particle and its neighbors, within a radius of influence, re. A weight function (kernel function), W(rij,re), determines the effect of a neighboring particle j (with position vector rj) on the target particle i (with position vector ri), where rij = |rj-ri| is distance between two particles. This study uses the third order polynomial kernel function, proposed by Shakibaeinia & Jin (2010) as: 
formula
(8)
A dimensionless parameter called particle number density, n, is used as a measure of density as: 
formula
(9)
here 〈 − 〉 is the weight averaging operator. The MPS approximation of the spatial derivatives (interpolation, gradient, divergence, and Laplacian) for the arbitrary scaler f and vector f on target particle i is given by (Shakibaeinia & Jin 2010): 
formula
(10)
where, n0 is the initial particle number density and λ is a coefficient defined by: 
formula
(11)
Figure 1

Schematic of MPS numerical approximation.

Figure 1

Schematic of MPS numerical approximation.

Pressure gradient

The pressure calculation method of this study is based on the weakly compressible MPS (WC-MPS) method (Shakibaeinia & Jin 2010), which assumes that the fluid is slightly compressible and uses an explicit EOS to calculate pressure from the density field. Tait's EOS (Shakibaeinia & Jin 2010) is used here, as: 
formula
(12)
where, c0 is an artificial (numerical) sound speed, selected to keep the compressibility (density change) lower than an acceptable value (typically less than 1%). The pressure gradient is then approximated using a modified form of MPS gradient formula as: 
formula
(13)
where

Rheological model implementation

As was mentioned, the model of this study considered the two-phase sediment water systems as a single multi-density multi-viscosity continuum. The fluid-type and sediment-type particles are used to discretize the continuum, where the fluid-type particles represent ambient water, and the sediment-type particles represent sediment grains and pore water between the grains (Figure 2). Each particle is assigned a density and viscosity. The water density is assigned to the fluid-type particles, and sediment bulk density, ρb, is assigned to the sediment-type particles.

Figure 2

Schematic of particle representation of sediment–water two-phase system.

Figure 2

Schematic of particle representation of sediment–water two-phase system.

The volume fraction (of the sediment phase) is explicitly calculated from the presence of the MPS particles. It is calculated by dividing the particle number density considering only the sediment phase (Ωs) particles, and the total particle number density as . This volume fraction is then used for determining the interface and calculation of the density field. The effective viscosity, η, of each particle is also needed, in order to calculate the divergence of shear stress tensor (; τ =2ηE). The effective viscosity of sediment-type particles is variable (a function of stress and strain rate) and obtained from the rheological model, while that of fluid-type particle is constant (η=ηf). From the μ(I) rheology formulation, discussed before, one can calculate the effective viscosity of a sediment-type particle i as: 
formula
(14)
The MPS approximation of two-dimensional strain-rate tensor, E, on particle i is given by: 
formula
(15)
The mechanical pressure between the sediment phase grains, ps, for incompressible dry granular flow is equal to the pressure p, given by the EOS. However, for submerged sediments (the case of this study), the pore pressure, ppore, must be subtracted. MPS method is known to be associated with some unphysical pressure fluctuations, which can affect yield behavior. To avoid this issue, some of the past studies (e.g., Shakibaeinia & Jin 2011; Manenti et al. 2011) have used a static pressure (assuming the vertical acceleration to be small negligible), as: 
formula
(16)
where, hi is vertical distance of particle i and the particle at sediment–water interface. However, this study uses a more generalized dynamic pressure, proposed by Jafari-Nodoushan et al. (2018), as: 
formula
(17)
This formula accounts for the vertical accelerations, and minimizes the effect of pressure fluctuations as it uses the smoothed value of thermodynamic pressure, 〈pi.
Now, the divergence of incompressible shear stress tensor can be calculated from the effective viscosity by: 
formula
(18)
As was mentioned, the MPS approximation of derivatives is based on interaction of neighboring particles. Considering the spatially variable viscosity field, when two neighboring particles i and j, with viscosities of ηi and ηj, are interacting, none of the viscosities can be used. Shakibaeinia & Jin (2011) proposed a harmonic mean of viscosities for such interaction, given by ηij = 2ηjηj/(ηi + ηj). Therefore, the MPS approximation of the divergence of shear stress tensor for particle i can be written as: 
formula
(19)

Boundary conditions

Free surface pressure is assigned to free surface boundaries. To locate the free surface particles, the value of particle number density is used. Particle number decreases while approaching the free surface. This is used as the criterion for recognition of free surface particles. A particle will be determined as a free-surface particle if its particle number density is less than β times of the initial average particle number density . In this paper, a β value of 0.94 (recommended by Shakibaeinia & Jin (2010)) is used.

For rigid boundaries (bed and walls), a few layers of so-called ghost particles are considered outside the boundary to assign the boundary values and to avoid unwanted density deficiency near the boundary. The pressure of ghost particles is extrapolated from the interior of the computational domain (to assure the boundary repulsive force). For the particles within the boundary layer, a logarithmic velocity profile is used. The details about the boundary treatment can be found in Jafari-Nodoushan et al. (2016).

Solution method

The time integration and solution method are based on a prediction correction algorithm, where each time step is divided into two sub-steps of prediction and correction. The velocity at a new time step (k + 1) is the sum of predicted velocity (u*) and corrected velocity (u′): 
formula
(20)
where 
formula
The relaxation factor α is chosen between 0 and 1 (here α= 0.5). In each time step the rheological model is applied to calculate the effective viscosity, which is then used to predict the velocity. Once the predicted velocity is calculated, particles are displaced, then the predicted value of the particle number density n* is calculated. The pressure is then computed from this particle number density. The pressure gradient term is then recalculated and used for calculation of corrected velocity and particle position. To account for the effect of high Reynolds flow, especially in water phase, a turbulence model is applied to calculate an eddy viscosity which is then added to the effective viscosity. Here, a large eddy simulation sub-particle scale (SPS) model (Shao & Gotoh 2005) is used. Details about the MPS implication of this model are available (Shakibaeinia & Jin 2010; Jafari-Nodoushan et al. 2016). This model has been used because it has been well-tested for particle methods. Time step size must satisfy the Courant–Friedrichs–Lewy (CFL) stability condition as well as stability of the diffusion equation as: 
formula
(21)

where dp is distance between the particles (size of the particles) and C is the Courant number and ηmax is maximum viscosity.

TEST CASES

In this section, first the model is validated for a benchmark, a viscoplastic Poiseuille flow, to evaluate the accuracy of the model in simulation of yielded and non-yielded regions in the flow field. Then, the developed model is validated and evaluated for simulation of sediment scouring in two common scouring situations: (1) sediment scouring under pipes and (2) sediment scouring induced by a water jet.

Viscoplastic Poiseuille flow

The geometry of the viscoplastic flow between two static parallel plates is shown in Figure 3(a). The flow is driven by the gravity, g. Note that flow occurs only if ρ|g| > y, where τy is a constant yield stress. Here, a dimensionless case with parameter values of h= 0.5, g =(0, −4/3), τy = μ1 × ps= 1, has been considered. The rheological model is assumed to be linear (I < <I0) to simplify the analytical solution as of a Bingham-plastic fluid and the material properties have been set to have equal to 0.5. The regularization exponent is set to m= 200. Domain is represented by particles with size of 0.0125. The steady-state analytical solution is taken from Chatzimina et al. (2005). In order to visualize the deformations, the particles are initially colored to form several parallel horizontal layers. Figure 3(b) shows a snapshot of the numerical results for the particle position at t= 2 (when time simulation reaches the steady condition). The particle layers show two distinguishable regions, yielded (near walls) and non-yielded (near center). The velocity numerical and analytical velocity profiles have also been compared showing a very good compatibility, with RSME of 0.008.

Figure 3

(a) The geometry of plane Poiseuille flow; (b) a snapshot of simulated final particle position (left), and comparison of numerical and analytical velocity profiles (right).

Figure 3

(a) The geometry of plane Poiseuille flow; (b) a snapshot of simulated final particle position (left), and comparison of numerical and analytical velocity profiles (right).

Sediment scour under circular pipe

In this section, the local scour beneath a submerged circular pipeline, based on the experiment of Mao (1986), is modeled. The sediment materials have a grain size of , bulk density of , and internal friction angle of degrees. The deposited layer has a thickness of 0.25 m, the water depth was 0.35 m and the pipeline diameter was 0.1 m. A two-dimensional computational domain with a length of 45 D and a height of 6 D is used in the simulation. The distance between the pipeline center and the inlet boundary is 15 D. The depth-averaged velocity at the inlet boundary is specified as 0.65 m/s. Figure 4 shows the general schematic of the forms of motion for the sand grain beneath the pipe.

Figure 4

The forms of motion for the sand grain beneath the pipe.

Figure 4

The forms of motion for the sand grain beneath the pipe.

The computational domain is represented with 28,508 particles (with initial distance of 0.01 m). The computational time on a machine with an Intel Core™ i5 (2.9 GHz) processor and 8GB of memory is about 5 hrs for 1 sec of simulations. Figure 5 shows results of numerical snapshots of local scour around a circular pipe problem at t= 1.5, 2.5, and 5 seconds. The high velocity flow under the pipe and flow circulation create a scour hole that grows as the time goes on, until an equilibrium is reached. The eroded materials are deposited further downstream and create a sediment ridge. Comparisons of the experimental and numerical snapshots show similar scouring and deposition features. A quantitative comparison of experimental and numerical sediment surface profiles (Figure 6) shows a good agreement, especially for the scouring zone. The dimensionless depth of scour hole (Ds = ds/D) is considered as the measure of the scouring.

Figure 5

Result of scour hole for MPS-WC model runs at t= 1.5, 2.5, and 5 seconds.

Figure 5

Result of scour hole for MPS-WC model runs at t= 1.5, 2.5, and 5 seconds.

Figure 6

Quantitative comparison of scour hole for MPS-WC model and experimental runs at t= 1.5, 2, and 5 seconds.

Figure 6

Quantitative comparison of scour hole for MPS-WC model and experimental runs at t= 1.5, 2, and 5 seconds.

The proposed numerical model has been able to successfully reproduce the shape and size of the scour hole. The numerical deposition zone, however, seems to be more shifted to upstream. This can be due to numerical assumption used for calculation of pore and mechanical pressure (which assumes a negligible exchange of water mass between ambient water and pore water). Figure 7 presents the time evolution of the dimensionless depth of scour hole. The experimental and numerical results are quite compatible. Figure 8 provides a sample snapshot of simulated, volume fraction, velocity magnitude, and velocity vectors at t= 1.5 sec. In this figure, the flow pattern, characterized by high-velocity flow under the pipe and circulation behind the pipe, can be seen. This erodes the sediment materials under the pipe and deposits them further downstream. The viscosity field differentiates the immobile particles (with very large viscosity) and particles in incipient motion (low-viscosity particles). These low-viscosity particles can be observed under the pipe, where the erosion happens. The present model has not used any empirical relation for calculation of erosion and deposition of sediments.

Figure 7

Comparison of numerical and experimental dimensionless scoured depth below the pipe center.

Figure 7

Comparison of numerical and experimental dimensionless scoured depth below the pipe center.

Figure 8

Simulated, volume fraction, velocity magnitude, and velocity vectors at t= 1.5 sec.

Figure 8

Simulated, volume fraction, velocity magnitude, and velocity vectors at t= 1.5 sec.

Jet-induced sediment scouring

Outflow from most hydraulic structures is in the form of a jet, which can be a source of severe sediment scouring at downstream channels. Here, the sediment scouring problem due to a non-submerged water jet (with various jet impact angles) in a rectangular plunge pool is investigated with the proposed mesh-free model of this study. The experimental data, required for comparison purposes, are taken from Pagliara et al. (2008). Figure 9 shows schematic geometry of the problem and Table 1 shows the geometrical parameters for two setups (A and B).

Figure 9

Schematic of jet scouring problem.

Figure 9

Schematic of jet scouring problem.

Table 1

Parameters of the experimental setup (Pagliara et al. 2008)

Setup h0/D Jet diameter D (m) Angle α (°) Discharge Q (L/s) 
5.1 0.03069 60 1.5 
7.2 0.03069 45 1.5 
Setup h0/D Jet diameter D (m) Angle α (°) Discharge Q (L/s) 
5.1 0.03069 60 1.5 
7.2 0.03069 45 1.5 

The water jet in these selected setups for modeling purposes has impact angles of 45 and 60° (with respect to a horizontal reference line). The sediment materials have a grain size of , bulk density of , and internal friction angle of . The depth (h0) and width of the channel and thickness of the sediment layer are 0.7 m, 0.5 m, and 0.3 m, respectively. Parameter is dimensionless time. Numerical modeling was performed for two setups, A and B, with jet impact angles of 45 and 60°. The solution domain was represented with particles with average distance (particle size) of 0.00115 m. Figure 10 compares the experimental and numerical snapshots (showing particle types) of jet scouring problem for 60° jet (setup A) at t= 3, 7, 15, and 47 seconds. The jet flow and the resulting circulations create a high-shear stress region (whose shear stress exceeds the yield stress), where the sediment material starts to be eroded. As a result, a scour hole forms and grows as the time goes on, until equilibrium is reached. The eroded materials are deposited further downstream under the gravity force and create a sediment ridge. Similarly, Figure 11 shows numerical snapshots of 45° jet (setup B) at t= 25, 50, 75, and 100 seconds. As expected, the 45° impact angle results in a less deep and more stretched scour hole compared to the 60° impact angle. Comparisons of the experimental and numerical snapshots show similar scouring and deposition features.

Figure 10

Snapshots of simulated and experimental (Pagliara et al. 2008) jet scouring problem (at t= 3, 7, 15, and 47 s for α = 60°).

Figure 10

Snapshots of simulated and experimental (Pagliara et al. 2008) jet scouring problem (at t= 3, 7, 15, and 47 s for α = 60°).

Figure 11

Snapshots of simulated jet scouring problem (at t= 5, 25, 50, and 100 s for α = 45°).

Figure 11

Snapshots of simulated jet scouring problem (at t= 5, 25, 50, and 100 s for α = 45°).

Figure 12 provides a quantitative comparison of experimental and numerical sediment surface profiles for the 60° jet (Case A). Note that for Case B, the experimental snapshots are not clear enough for accurate extraction of the profiles. Figure 12 shows a good compatibility of the profiles, especially for the scour hole. The proposed numerical model has been able to successfully reproduce the shape and size of the scour hole.

Figure 12

Comparison of experimental and numerical bed profiles at t = 7, 15, and 47 seconds for α = 60°.

Figure 12

Comparison of experimental and numerical bed profiles at t = 7, 15, and 47 seconds for α = 60°.

Two parameters, including the dimensionless depth of scour hole (Ds = ds/D) and dimensionless height of deposition ridge (Dd = dd/D), are considered as the measures of the scouring and deposition features. Figure 13 presents the time evolution of Ds and Dd for jet angles (α= 60 and 45°). The experimental and numerical results are quite compatible (especially for the scour hole depth) for both jet angles. The numerical model slightly underestimates the deposition ridge height as it stretches the deposition region to the downstream. The root mean square error provided in Figure 14 also shows a better prediction of the scour hole depth compared to the deposition ridge height. Accuracy of the results for the sediment deposition can be due to the numerical assumption used for calculation of pore and mechanical pressure (which is not suitable for suspended materials), as well the two-dimensional simulation of a three-dimensional case.

Figure 13

Comparison of numerical and experimental dimensionless scour depth (Ds) and dimensionless depth at the outlet of scour hole (Dd).

Figure 13

Comparison of numerical and experimental dimensionless scour depth (Ds) and dimensionless depth at the outlet of scour hole (Dd).

Figure 14

Graph of measured and simulated dimensionless depth of scour hole parameters.

Figure 14

Graph of measured and simulated dimensionless depth of scour hole parameters.

Figure 15 provides a sample snapshot of simulated viscosity, velocity magnitude/vector, and volume fraction at t= 5 s. The flow circulation pattern has also been highlighted in this figure. As the jet flow hits the sediment surface, it is divided into two branches: (1) a spoon-like counterclockwise strong flow at the direction of the jet and (2) a weak backward (clockwise) flow in the opposite direction of the jet. These flow currents erode the sediment material, creating a scour hole, and deposit them downstream and upstream, respectively. As time goes on, the scour hole grows. Where the shear stress due to these flow currents exceeds the yield stress (i.e., low viscosity regions) the sediment materials start to move and separate from the bed. The material is carried by the flow current and deposited further downstream under their weight. This creates an active (yielded) sediment layer on the sediment surface. The sediment layer under this active layer remains relatively stagnant (un-yielded). This is because of the relatively low shear stress, high yield stress, and effective viscosity (due to the weight of the above layers).

Figure 15

Simulation of free surface profile, values of velocity, velocity vectors, and volumetric fraction of sediments at t= 5 seconds and α = 45°.

Figure 15

Simulation of free surface profile, values of velocity, velocity vectors, and volumetric fraction of sediments at t= 5 seconds and α = 45°.

CONCLUSIONS

Multiphase mesh-free particle modeling of local sediment scouring with μ(I) rheology was developed for continuum modeling and evaluated for simulation of sediment scouring. The sediment material is considered as a non-Newtonian viscoplastic fluid whose behavior is predicted using a regularized μ(I) rheological model in combination with pressure-dependent yield criteria. The model is first validated for a benchmark problem of viscoplastic Poiseuille flow. Then, the developed model is validated and evaluated for two classical sediment scouring cases, namely, scouring under submerged pipes, and scouring resulting from water jets with different angles. The proposed numerical model has been able to successfully reproduce the shape and size of the scour hole. The numerical deposition zone, however, seems to be more shifted to upstream. This can be due to the numerical assumption used for the calculation of pore and mechanical pressure. The experimental and numerical results for the jet are quite compatible (especially, for the scour hole depth) for both jet angles. The numerical model slightly underestimates the deposition ridge height as it stretches the deposition region to the downstream. The results show how the highly erosive flow and circulations can create a low-viscosity region at the surface of the sediment continuum, where the shear stress can exceed the yield stress and the sediment can be eroded to form the scour hole. The results of this study not only demonstrated the capability of mesh-free particle methods for sediment scouring modeling, but also provided a better understanding of flow and sediment scouring in the test cases of this study.

REFERENCES

REFERENCES
Ataie-Ashtiani
B.
&
Farhadi
L.
2006
A stable moving-particle semi-implicit method for free surface flows
.
Fluid Dynamics Research
38
,
241
256
.
Azamathulla
H. M.
&
Ghani
A. A.
2010
Genetic programming to predict river pipeline scour
.
Journal of Pipeline Systems Engineering and Practice
1
(
3
),
127
132
.
Barati
R.
,
Neyshabouri
S. A. A. S.
&
Ahmadi
G.
2018
Issues in Eulerian–Lagrangian modeling of sediment transport under saltation regime
.
International Journal of Sediment Research
33
(
4
),
441
461
.
Calantoni
J.
,
Holland
K. T.
&
Drake
T. G.
2004
Modelling sheet-flow sediment transport in wave-bottom boundary layers using discrete-element modeling
.
Philosophical Transactions-Royal Society of London, Series A, Mathematical, Physical and Engineering Sciences
362
,
1987
2002
.
Capone
T.
,
Panizzo
A.
&
Monaghan
J. J.
2010
SPH modelling of water waves generated by submarine landslides
.
Journal of Hydraulic Research
48
(
S1
),
80
84
.
Chatzimina
M.
,
Georgiou
G. C.
,
Argyropaidas
I.
,
Mitsoulis
E.
&
Huilgol
R. R.
2005
Cessation of Couette and Poiseuille flows of a Bingham plastic and finite stopping times
.
Journal of Non-Newtonian Fluid Mechanics
129
(
3
),
117
127
.
Chen
C. L.
&
Ling
C. H.
1996
Granular-flow rheology: role of shear-rate number in transition regime
.
Journal of Engineering Mechanics
122
(
5
),
469
480
.
Fadafan
M. A.
&
Kermani
M. R. H.
2017
Moving particle semi-implicit method with improved pressures stability properties
.
Journal of Hydroinformatic
20
(
6
),
1268
1285
.
Falconer
R. A.
&
Owens
P. H.
1990
Numerical modelling of suspended sediment fluxes in estuarine waters
.
Estuarine, Coastal and Shelf Science
31
(
6
),
745
762
.
Forterre
Y.
&
Pouliquen
O.
2008
Flows of dense granular media
.
Annual Review of Fluid Mechanics
40
,
1
24
.
Gingold
R. A.
&
Monaghan
J. J.
1977
Smoothed particle hydrodynamics: theory and application to non-spherical stars
.
Monthly Notices on the Royal Astronomical Society
181
,
375
389
.
Gotoh
H.
&
Sakai
T.
2006
Key issues in the particle method for computation of wave breaking
.
Coastal Engineering Journal
53
(
2–3
),
171
179
.
Guan
D.
,
Melville
B.
&
Friedrich
H.
2016
Local scour at submerged weirs in sand-bed channels
.
Journal of Hydraulic Research
54
(
2
),
172
184
.
Hsu
T. J.
,
Jenkins
J. T.
&
Liu
P. L. F.
2003
On two-phase sediment transport: dilute flow
.
Journal of Geophysical Research (Oceans)
108
,
14
.
Hu
X. Y.
&
Adams
N. A.
2007
An incompressible multi-phase SPH method
.
Journal of Computational Physics
227
(
1
),
264
278
.
Jafari-Nodoushan
E.
,
Hosseini
K.
,
Shakibaeinia
A.
&
Mousavi
S. F.
2016
Meshless particle modelling of free surface flow over spillways
.
Journal of Hydroinformatics
18
(
2
),
354
370
.
Jafari Nodoushan
E.
,
Shakibaeinia
A.
&
Hosseini
K.
2018
A multiphase meshfree particle method for continuum-based modeling of dry and submerged granular flows
.
Powder Technology
335
,
258
274
.
Jop
P.
,
Forterre
Y.
&
Pouliquen
O.
2006
A constitutive law for dense granular flows
.
Nature
441
(
7094
),
727
730
.
Kashyap
S.
,
Dibike
Y.
,
Shakibaeinia
A.
,
Prowse
T.
&
Droppo
I.
2017
Two-dimensional numerical modelling of sediment and chemical constituent transport within the lower reaches of the Athabasca River
.
Environmental Science and Pollution Research
24
(
3
),
2286
2303
.
Khanpour
M.
,
Zarrati
A. R.
,
Kolahdoozan
M.
,
Shakibaeinia
A.
&
Amirshahi
S. M.
2016
Mesh-free SPH modeling of sediment scouring and flushing
.
Computers & Fluids
129
,
67
78
.
Koshizuka
S.
&
Oka
Y.
1996
Moving particle semi-implicit method for fragmentation of incompressible fluid
.
Nuclear Science and Engineering
123
(
3
),
421
434
.
Larese
A.
,
Rossi
R.
,
Oñate
E.
,
Toledo
M. Á.
,
Morán
R.
&
Campos
H.
2013
Numerical and experimental study of overtopping and failure of rockfill dams
.
International Journal of Geomechanics
15
(
4
),
04014060
. doi: 10.1061/(ASCE)GM.1943-5622.0000345
Manenti
S.
,
Sibilla
S.
,
Gallati
M.
,
Agate
G.
&
Guandalini
R.
2011
SPH simulation of sediment flushing induced by a rapid water flow
.
Journal of Hydraulic Engineering
138
(
3
),
272
284
.
Mao
Y.
1986
The Interaction Between a Pipeline and an Erodible bed
.
Ph.D. thesis
,
Technical University of Denmark
,
Lyngby
,
Denmark
.
Nabian
M. A.
&
Farhadi
L.
2016
Multiphase mesh-free particle method for simulating granular flows and sediment transport
.
Journal of Hydraulic Engineering
143
(
4
),
04016102
.
Doi: 10.1061/(ASCE)HY.1943-7900.0001275
Pagliara
S.
,
Hager
W. H.
&
Unger
J.
2008
Temporal evolution of plunge pool scour
.
Journal of Hydraulic Engineering
134
(
11
),
1630
1638
.
Papanastasiou
T. C.
1987
Flows of materials with yield
.
Journal of Rheology
31
(
5
),
385
404
.
Park
J. H.
,
Sok
C.
,
Park
C. K.
&
Do Kim
Y.
2016
A study on the effects of debris accumulation at sacrificial piles on bridge pier Scour: I. Experimental results
.
KSCE Journal of Civil Engineering
20
(
4
),
1546
1551
.
Roushangar
K.
,
Mirheidarian
S.
&
Shiri
J.
2017
Modeling local pier scour with bed effect implications: heuristic vs. empirical strategies
.
ISH Journal of Hydraulic Engineering
23
(
1
),
13
22
.
Rycroft
C. H.
,
Wong
Y. L.
&
Bazant
M. Z.
2010
Fast spot-based multiscale simulations of granular drainage
.
Powder Technology
200
(
1–2
),
1
11
.
Shakibaeinia
A.
&
Jin
Y. C.
2010
A weakly compressible MPS method for simulation of open-boundary free-surface flow
.
International Journal for Numerical Methods in Fluids
63
(
10
),
1208
1232
.
Shakibaeinia
A.
&
Jin
Y. C.
2011
A mesh-free particle model for simulation of mobile-bed dam break
.
Advances in Water Resources
34
,
794
807
.
Shakibaeinia
A.
&
Jin
Y. C.
2012a
MPS mesh-free particle method for multiphase flows
.
Computer Methods in Applied Mechanics and Engineering
229–232
,
13
26
.
Shakibaeinia
A.
&
Jin
Y. C.
2012b
Lagrangian multiphase modeling of sand discharge into still water
.
Advances in Water Resources
48
,
55
67
.
Shao
S.
&
Gotoh
H.
2005
Turbulence particle models for tracking free surfaces
.
Journal of Hydraulic Research
43
(
3
),
276
289
.
Tajnesaie
M.
,
Shakibaeinia
A.
&
Hosseini
K.
2018
Meshfree particle numerical modelling of sub-aerial and submerged landslides
.
Computers & Fluids
172
,
109
121
.
Tao
J.
&
Li
J.
2018
CFD-DEM Two-Way Coupled Numerical Simulation of Bridge Local Scour Behavior Under Clear-Water Conditions
.
Transportation Research Board 97th Annual Meeting (No. 18-05939)
,
Washington, DC
,
USA
.
van Rijn
L. C.
1984
Sediment transport, part I: bed load transport
.
Journal of Hydraulic Engineering
110
(
10
),
1431
1456
.