Multiphase mesh-free particle modeling of local sediment scouring with μ(I) rheology

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) developed 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 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.

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.

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).

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

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.

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| > hτ_y, where τ_y is a constant yield stress. Here, a dimensionless case with parameter values of h= 0.5, g =(0, −4/3), τ_y = μ₁ × ps= 1, has been considered. The rheological model is assumed to be linear (I < <I₀) 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.

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.

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 hours for 1 second 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 (D_s = d_s/D) is considered as the measure of the scouring.

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 seconds. 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.

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).

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 (h₀) 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 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.

Two parameters, including the dimensionless depth of scour hole (D_s = d_s/D) and dimensionless height of deposition ridge (D_d = d_d/D), are considered as the measures of the scouring and deposition features. Figure 13 presents the time evolution of D_s and D_d 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 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).

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.