Mesh-free SPH modeling of sediment scouring and flushing

Computers and Fluids 129 (2016) 67–78

Contents lists available at ScienceDirect

Computers and Fluids journal homepage: www.elsevier.com/locate/compfluid

Mesh-free SPH modeling of sediment scouring and flushing M. Khanpour a, A.R. Zarrati a, M. Kolahdoozan a, A. Shakibaeinia b,∗, S.M. Amirshahi a a b

Department of Civil and Environmental Engineering, Amirkabir University of Technology, Tehran, Iran Water & Climate Impact Research Centre, University of Victoria, Victoria, Canada

a r t i c l e

i n f o

Article history: Received 24 February 2015 Revised 27 September 2015 Accepted 6 February 2016 Available online 15 February 2016 Keywords: Smoothed Particle hydrodynamics (SPH) Continuum-based sediment transport Sediment scouring Sediment flushing

a b s t r a c t A two-phase mesh-free Lagrangian model, based on the Smoothed Particle Hydrodynamics (SPH) method, is developed and evaluated for the continuum-based study of sediment-water system (a dense submerged granular flow system). The case studies are (a) sediment scouring downstream of a wall-jet, and (b) reservoir sediment flushing, both of great interest of engineers and researchers. The mesh-free Lagrangian nature of the SPH method makes the model capable of dealing with the large interfacial deformations involved in these complex water-sediment flow systems. Water and sediment materials are both treated as continuum and the governing equation of their motion are solved on a Lagrangian frame using SPH formulation. A pressure-dependent viscoplastic rheological model is used to describe the behavior of the sediment phase. A combination of Mohr–Coulomb and Shields yielding criteria (in place of commonly used Mohr–Coulomb yield criteria) is employed to determine the critical stress above which the sediment material behaves as a viscous fluid. Effect of this combination on the model results is investigated. To validate the developed numerical model, this study also provides some original experimental measurements for the test cases of this study. The results show the good compatibility of the numerical and experimental water and bed profiles. The sediment feature such as scouring and deposition areas (for wall-jet scouring case) and the sediment failure slope (for the flushing test case) are reproduced by the developed SPH model, with a good degree of accuracy. The numerical results as well as experimental measurements of this study not only provide an insight into better understanding of the complicated mechanism of sediment transport induced by rapid flow discharge, but also provide a reliable numerical tool and benchmark for simulation of such problems. © 2016 Elsevier Ltd. All rights reserved.

1. Introduction The reliable estimation of sediment transport is one of the most importance tasks in open-channel hydraulic and related areas. The prediction of sediment transport especially for case of turbulent highly erosive and transient flows can be very challenging. Until now, prediction of flow and sediment transport in open-channels has been mostly based on Eulerian mesh-based numerical methods such as finite element and finite volume (e.g. in [1–3]). Eulerian methods solve flow equation (conservation of mass and momentum) over a mesh system. These methods take care of sediment transport typically by solving the transport equations (for suspended sediments) and different empirical equations for bedload materials. Every empirical equation has its own conditions and limitations. Thus, these methods can be accompanied by uncertainties in determining parameters such as critical shear stress [4–6], settling velocity [7,8] reference concentration [9–11] bed

∗

Corresponding author. Tel.: +1250 363 8963. E-mail address: [email protected], [email protected] (A. Shakibaeinia).

http://dx.doi.org/10.1016/j.compfluid.2016.02.005 0045-7930/© 2016 Elsevier Ltd. All rights reserved.

layer thickness [12], bed load [9,13], and the non-equilibrium adaptation length and coefficients [14–17]. Furthermore, mesh-based methods have problem dealing with the large interfacial deformations (both at the free surface and water-bed interface) which is a common feature of highly erosive turbulent flows (Shakibaeinia and Jin [18,20]. The most recent generation of numerical methods, the meshfree Lagrangian methods (also called particle methods) have become a serious alternative for the traditional mesh-based methods especially for the applications with large interfacial/boundary deformations and fragmentations. Smoothed Particle Hydrodynamics (SPH), developed by Gingold and Monaghan [21] and Lucy [22], and Moving Particle Semi-implicit Method (MPS) [23] are two of the widely-used particle methods. SPH is being increasingly used to model various types of fluid flows including, viscous flows [24,25], free surface flows [26,27], flow through porous media [28,29], heat conduction [30], and explosion [31]. SPH has also been developed and applied for multiphase flows (e.g. in [32–35]). Monaghan and Kocharyan [32] formulated the two-phase flow of a dusty gas produced during the eruption of a volcano using SPH. . Colagrossi and Landrini [33] studied

68

M. Khanpour et al. / Computers and Fluids 129 (2016) 67–78

two-dimensional interfacial flows to simulate sediment transport; in which, water and sediments were treated a two phase system of Newtonian and non-Newtonian. The similar approach was also used in [[36–38], and [39]]. Ataie-Ashtiani et.al. [40] also simulated land-slide by this method. Zanganeh et.al. [41] simulated scour under the pipeline by this technique in combination with DEM model. Shakibaeinia and Jin ([18–20]) proposed a MPS based model for multiphase granular flows, and evaluated various rheological models for simulations of sediment transport. In all of these studies were based on continuum description of sediment phase. A constitutive model (typically Bingham Plastic) has been used to describe the behavior of sediment particles whose shear stress exceeds a yielding criteria (e.g., Mohr–Coulomb). In this study a SPH model is developed for the multiphase systems of water and sediments in the case of sediment scouring and flushing. The model treats the sediment phase as a continuum (i.e. a visco-plastic fluid). The commonly used Mohr–Coulomb yield criterion is replaced with a combination of Mohr–Coulomb and Shields in order to (1) achieve a more reliable erosion pattern; (2) relate the rheological behavior of sediment materials to the commonly measured sediment properties (e.g. sediment grain size). The algorithms for calculation of effective pressure and nearbed velocity profiles are also developed. Effect of these modifications on the predicted sediment profiles is evaluated for the test cases of this study. Two sets of original experiments on a sediment scouring (resulted from a wall-jet) and reservoir flushing test cases are conducted and their results are used to validate the developed numerical SPH model. This study is expected to result in the development of an accurate mesh-free Lagrangian numerical tool for modeling sediment transport. It will also provide a better insight into mechanisms involved in transport of sediments for cases wall-jet scouring and reservoir flushing using a combination of numerical and experimental analysis. 2. The numerical model To model the multiphase system, we consider the isothermal Navier–Stokes equations on a Lagrangian frame given by:

Dρ + ρ (∇ · u ) = 0 Dt Du ρ = −∇ p + ∇ (μ∇ · u ) + F Dt

(1)

where ρ is the density, u is the flow velocity, F represents the body forces (e.g. gravity g), μ is the dynamic viscosity, and p is the thermodynamic pressure. There is no convective acceleration term in the Lagrangian frame momentum equation, and the movement of particles is simply given by Dr/Dt = u, with r being the position vector. Here the Tait’s equation of state (Batchelor, 1967) is used to explicitly calculate the thermodynamic pressure.

p=B

 ρ γ ρ0

−1



(2)

The typical value of γ = 7 is used in this relation. The B parameters in this equation are selected using scale analysis [[42,24,25], and [43]] so that the density variation is less than a certain value. 2.1. SPH method The smoothed-particle hydrodynamics (SPH) represents the material domain by a set of discrete elements with no connections (referred to as particles), over which the flow equations are integrated using a kernel (smoothing) function. The arbitrary quantity A(r) on the particle i are approximated by summing up the values

Fig. 1. Sediment particles and fluid particles and interface domain in two phase flow.

on the neighboring particles j given by (Liu [44]):



A(r )i =

mj

j

Aj

ρj

W ( |ri − r j |, h )

(3)

where the operator   is kernel smoothing operator, r is the position vector, m is mass of particle, W is kernel (or smoothing) function, and h is smoothing length. Throughout the present simulations, the commonly used cubic spline kernel function is used [44]. The SPH spatial derivative of a quantity can be obtained by:

∇ A(r )i =



mj

j

Aj

ρj

∇ W ( |ri − r j |, h )

(4)

In the original SPH method, the density approximation of a particle is expressed by summing the mass contributions from the neighboring particles densities ([24]). Employing this approach can lead to some oscillation near the water-sediment interface due to the density gradient (Fig. 1). To overcome this problem, here the densities (in both water and sediment phases) are calculated by (Hu and Adams, [34]):

ρi = m i



W ( |ri − r j |, h )

(5)

j

In comparison to the commonly used formulation which relates the density of the particle of interest to the mass of its neighboring particles ( j mj Wij ), in Eq. (5), the particle density depends only on the local distribution of the particles and therefore performs better near the interfaces [34]. The system treated as multi-density multi-viscosity system. The SPH discretization of the momentum equation, used in this study for both phases, is given by:

 Du  Dt

i

=

F

ρi ×

−





mj

j

ui j .ri j r2i j + 0.01h

  m j 4 μi μ j + ∇iWi j + ρi ρ j μ i + μ j ρi2 ρ 2j j pi

∇iWi j

pj

(6)

where uij = uj -ui and rij = rj -ri are the relative velocity and position vector, respectively. The first and second summations (in the left hand side of the equation) estimate the pressure gradient and viscous term, respectively. The viscous term uses an interaction viscosity equal to the harmonic mean of μi and μj viscosities as suggested by [19,34,45] for multiphase flow problems. The particles are moved based on the calculated velocity fields. The XSPH model introduced by Monaghan [42] is used to smooth the velocity field to increase the stability and reduce the un-physical fluctuations. In XSPH a velocity, corrected using the contribution from neighboring particles, is used to move particles in each time step as:

 mj d xi = ui − ε u W dt ρj ij ij

(7)

M. Khanpour et al. / Computers and Fluids 129 (2016) 67–78

69

Fig. 2. Classification of the sediment particle based on local distribution in the numerical model.

where ε is a constant in the range of 0<ε <1. Here in this study, an improved XSPH model (IXSPH) [27] is applied, where, the corrected velocity (with relatively smaller ε of ∼0.02) is used both in momentum equation and for moving the particles [27,46]. This modification is only used for water particles, because, the highly viscous nature of the granular phase automatically damps the velocity fluctuation, resulting a smooth velocity field. To simulate the solid boundaries, a layer of boundary particles are used [42]. To resolve the near-boundary density deficiency, some layers of virtual particles are placed at the boundaries exterior. The pressures for both phases are obtained from the Tait’s equation of state (2). Unlike for the water phase particles, the dynamic viscosity for the sediment phase particles is unknown, and need to be determined from a rheological model. 2.2. Sediment phase rheology The viscosity of sediment particles was calculated based on Bingham Plastic model, which is a simple yet effective and widely used model for granular flows [18,19,37,39,41, and 43]. Bingham Plastic behaves as a rigid body at stresses lower than a yield stress and flows as a viscous fluid when stresses exceed this yield stress. The constitutive relation of the Bingham plastic model is given by:

τ = 2 μe f f E ; μe f f =



τy / 4IIE + 2μ0 ∞

| τ | > τy | τ | ≤ τy

(8)

in which, μeff is effective viscosity, μ0 is Bingham viscosity, E = 0.5 (∇ u + (∇ u)T ) is strain rate tensor, |A| = IIA = (0.5 A:A) is the magnitude of the arbitrary tensor A (IIA being second invariant of the tensor), and τ y is the yield stress. Eq. (8) is singular and discontinues when |τ | approaches τ y . To avoid this, a regularized bi-viscous version of Bingham plastic model, in which the effective viscosity is limited to a sufficiently large maximum viscosity μmax , is used. The yield stress τ y is typically defined by Mohr–Coulomb yielding criteria for granular material as [18,19]:

τy = c cos ϕ + p sin ϕ

(9)

where c is the cohesion coefficient, ϕ is the internal friction angle and p’ is the effective or mechanical pressure (normal stress between sediment grains). Note that μ0, ϕ and c are the material properties (a function of parameters such as density, and grain size/shape) and are usually measured experimentally (i.e. using rheometres). In the present study the sediment material is assumed cohesionless (c = 0). The second invariants of stress tensor

and strain rate tensor are given by:

IIτ =



 1 2 2 2 + τ23 + τ31  2 τi j τi j − (τkk )2  = τ12 − (τ11 τ22 + τ11 τ33 + τ22 τ33 )

IIE =



 1  ei j ei j − (ekk )2  = e212 + e223 + e231 2 − (e11 e22 + e11 e33 + e22 e33 )

(10)

2.2.1. Implementation of rheological model The components of strain rate tensor, E, are calculated by: n m j β ∂ Wi j ∂ vi β  = v ∂ xα ρ j ji ∂ xi α j=1

(11)

Effective pressure (or mechanical pressure) is defined as the normal stress between sediment grains. For dry sediment (negligible pore pressure) it related to the thermodynamic pressure as p’ = p-λ∇ .u, where λ is the bulk viscosity. For cases, where the vertical acceleration of the granular materials is small enough, a static pressure (equal to the lithostatic/overburden pressure for dry sediments) can be used. For the submerged sediment the pore water pressure must be subtracted; therefore, the effective pressure of each sediment particle is given by the buoyant weight of the overlaying sediments.:

p i = (ρs − ρw )gZi g g

(12)

where Zi is the distance (in the direction of gravity vector) between sediment particles and the water-sediment interface, ρ w and ρ s are water and sediment (bulk) densities, respectively. To obtain Zi g at first the interface must be identified in each time step. The sediment particles can be classified into three groups based on their local distribution as: 1- suspended particles 2- bed sediment (or hidden) particles 3- interface particles (Fig. 2). Since suspended particles are surrounded by water particles (not in direct contact with other sediment particles) their yield stress will be zero (as φ = 0) and they flows as a viscous material [18] and [19]). To search the interfacial particles, the solution domain is divided to the equal-distance regions (columns) perpendicular to the direction of the gravity vector. The distance (≥ particle size) is selected based on the required accuracy (resolution) for interface detection. A representative interfacial particle is determined for each of these regions using several criteria. The interfacial particle has the highest elevation among the sediment type particles whose volume fraction is greater than a ξ value (its density is greater

70

M. Khanpour et al. / Computers and Fluids 129 (2016) 67–78

Fig. 3. Identifying the nearest water particles (j) to the interface particle (i).

than ξ times of densities of sediment+water), otherwise it belongs to the suspended sediment or hidden bed sediment particles.:

ρti ≤ ξ (ρ0s + ρ0w ) → suspended particles ρti > ξ (ρ0s + ρ0w ) → interface particles < ρt i >=

n 

m jWi j

; (13)

j=1

where ρ 0s is reference density of sediment particles. The ξ parameter defines whether a sediment particle is suspended or it is in contact with the other sediment particles. Theoretically, for uniformly distributed sediment-water mixture, a sediment particle whose smoothed (interpolated) density is less than half of sediment+water initial density (ξ = 0.5), has no contact with the other sediment particles and therefore can be considered as suspended particle. In practice, the numerical results are not very sensitive to this choice because the interaction viscosity of the sediment particles in Eq. (6) is defined by the harmonic mean of viscosity of that particle and its neighboring particle. Therefore, the value of this interaction viscosity will be very small when the sediment particle is surrounded by a high number of water particles so the sediment particle will yield anyway, regardless of the choice of ξ in such cases. To also rule out the overhanging sediment materials another criteria is considered, in which there must be no particle with a volume fraction close to zero under the detected interfacial particle. Note that this only happens if a group of water particles are located under the initially detected interface; in other hand the overhanging criteria does not imply for the single scattered water particles that might be trapped between the sediment particles). An efficient search algorithm is implied, which only requires one time search of the particle list. The particle list is searched one by one and if a particle satisfies the interfacial particle criteria it is considered as an interfacial candidate for the region (column) in which it is located. That particle becomes the interfacial representative of that particular column, unless a more appropriate candidate (e.g., the one with higher elevation) is found as the search continues. By the end of search each column will have an interfacial representative. After the interface particles are recognized, for each region, the Zg value can be calculated. To find Zg for a particle i, first the region to which the particle i belongs is determined, then the vertical distance of i to the interfacial representative of that region is calculated. Effective pressure is then calculated using Eq. (12). The sediment particles which are further away from the water interface (deeper at the bed) have a lower strain rate, higher effective pressure, yield stress and effective viscosity. The effective pressure and Mohr–Coulomb yield stress (and consequently the effective viscos-

Fig. 4. correction of the shear stress based on the bed slope.

ity) for approach zero for the near-interface sediment particles. Therefore, these particles cannot resist against the hydrodynamic forces. To overcome this problem an additional Shields’ yielding criteria is employed. The Shields’ critical shear stress is compared against a simplified bed shear stress, calculated using the boundary layer theory, to determine if an interfacial sediment particle can move. This also helps to related the rheological behavior of the sediment material to the commonly measure sediment properties such as grain size. To implement the boundary layer theory, the distance between the nearest water particle, j, to a sediment interface particle, i, (y) needs to be determined (Fig. 3). Knowing the stream-wise velocity of water particle j, streamwise velocity u, shear velocity, u∗ , and the shear stress, τ , on the sediment particles can be calculated as:

u

u∗ u u∗

= κ1 ln 1

= κ ln

30y ks

9u∗ y υ

Re∗ > 70 (Rough bed ) Re∗ < 5 (Smoothed bed )

τ = ρ u∗ 2

(14)

(15)

where κ = 0.4, ks is equivalent roughness of bed material (∼2d, d being the bed grain size) and Re∗ = ρ u∗ d/μ is the shear Reynolds number and μ is water viscosity. Note that the result of the numerical model shows that the bed condition was hydraulically rough. Shields parameter,θ , is given by (Van Rijn, [47]:



θ=

0.010595 ln(Re∗ ) + 0.110476/Re∗ + 0.0027197 0.068

for Re∗ ≤ 500 for Re∗ > 500

(16)

The Shields’ yield shear stress for sediment particles, τ s ,is then calculated by:

τs = θ (γs − γw )d

(17)

When the particle shear stress exceeds the Shields’ critical shear stress (τ >τ s ) the interfacial sediment particles are allowed to move as a viscous fluid. The critical shear stress derived by Eq. (17) is for a horizontal bed and need to be corrected by the bed slope, β , (Fig. 4):

M. Khanpour et al. / Computers and Fluids 129 (2016) 67–78

71

Fig. 5. Experimental setup for (a) wall-jet scouring, and (b) sediment flushing test case. Table 1 Material properties for the reservoir flushing and wall-jet scouring test cases. Sediment phase Sediment material Natural sand

τs = τs,0

Water phase Grain size d50 (mm) 0.85

Grain density ρ g (kg/m3 ) 2650

sin(φ ± β ) sin(φ )

Bulk density ρ s (kg/m3 ) 1750

(18)

The slope of interface at any point is locally calculated considering the position of each interface particle, i, in respect to its two nearest interface neighboring particles, j and k.

αik = tan

−1

y − y  i k xk − xi



αi j = tan

−1

y j − yi x j − xi



αi =

αi j + αik 2

Friction angle φ (°) 35

Density ρ w (kg/m3 ) 10 0 0

viscosity μ (pa.s) 0.001

15 cm (Fig. 5b). The gate is instantly opened and sand material is washed downstream from the gate. For both cases, a deposit of non-cohesive uniform sand (d50 = 0.85 mm), with the bulk density of ρ s = 1750 kg/m3 , and friction angle of ∼35° is used. Table 1. summarizes the material properties for both test cases. 3. Results and discussion

(19)

2.3. Experimental setup To verify the numerical model of this study, we conduct two sets of experiments on the small scale (1) sediment scouring (resulted from a wall-jet), and (2) reservoir sediment flushing. The experiments are done in a glass-walled horizontal flume of 0.74 m wide and 0.20 m deep (Fig. 5). A rigid plate (with a removable gate at its bottom portion) is installed in the flume to create an artificial reservoir in one side. The gate can rotate around its hinge and open instantly creating a bottom outlet for the reservoir. For wall jet scouring test, a 5 cm sand layer is laid downstream of the gate and leveled horizontally (Fig. 5a). The reservoir is filled with water with height of h = 15 cm. The gate is opened instantly and the scouring process is recorded. An image processing technique is then used to determine the water and bed elevation change in time. For the reservoir flushing test, one layer of sandy soil with a height of 5 cm is placed on the reservoir bed and leveled horizontally. The reservoir is filled with water with height of h =

3.1. Wall-jet scouring test case Water particles, sediment particles, and boundary/ ghost particles are placed with same space of 1.7 mm, (about two times of the grain size) resulting in total 46,160 particles. Needless to say, in the particle methods, particles are the points carrying the physical properties of the continuum, therefore particle spacing (or size) is should not be necessarily equal to the sediment grain size. Fig. 6. shows the configuration of the SPH model setup. After the simulation started the gate will be kept closed until the initial pressure fluctuations are damped and a relatively hydrostatic pressure distribution is established in the reservoir. At next stage, the boundary and ghost particles at the gate position are removed instantaneously resulting in formation of a wall jet and transport of sediment particles. The computation time of the model for this test case was ∼7 sec per time step for on a 2.4 GHz Intel i7 processor. Fig. 7. compares some snapshots of the numerical and experimental results. Fig. 8. compares the evolution of the experimental and numerical water and bed surface profiles. Once the gate is opened, a scour hole starts to form at the beginning of the channel as the result of the wall jet. Washed sediment particles are

72

M. Khanpour et al. / Computers and Fluids 129 (2016) 67–78

Fig. 6. SPH particle configuration for the wall-jet scouring test case.

Fig. 7. Sample snapshots of numerical and experimental results for the wall-jet scouring test case.

deposited further downstream in the form of a ridge, in the area where the bed shear stress is low enough. (t = 0.9 s). This deposition ridge is then transported downstream gradually. As the time goes on, the water level in the reservoir and therefore the water velocity decrease, thus the scouring process stops (at ∼2 s). The curvature of the free surface and bed surface are in the same phase both in experiments and the numerical model. At the end of ex-

periments a returning wave forms over the sediment ridge which is also correctly reproduced by the numerical model. By decreasing the water level in the reservoir the velocity of the water particles are reduced. Thus the ability of these particles to pass the water surface hump (a wave) is diminished resulting in an increase in the height of the water surface hump (t = 1.9 s). Subsequently the slope of the upstream side of the hump increases

M. Khanpour et al. / Computers and Fluids 129 (2016) 67–78

73

Fig. 8. Evolution of experimental and numerical water and bed profile features for wall-jet scouring test case.

and the free surface profile becomes unstable (t = 2.1 s). On the other hands as the times goes on, the water level on the hump became higher than the water level in the reservoir. This creates a surface wave propagating upstream toward the reservoir (t = 2.1 s to 3.1 s). As it is illustrated the present numerical method can successfully reproduces the complex water and bed surface of this test case proving the capabilities of the present numerical model. The numerical and experimental profiles show a good compatibility. The free surface and bed profile features are similar. There are however some small discrepancies. For instance, the water level at t = 1.4 s reached under the gate level while the free surface in the numerical method is still slightly higher than the gate. Part of the discrepancy can be due to the mechanism of opening of the gate (numerical model has an instant gate removal while at the experiment the gate is removed within a time frame of ∼0.2 s. The effect of the proposed yield stress calculation method, which uses a combination of Mohr-Coulomb τ y and Shields τ s

stresses, in comparison with the commonly used Mohr–Coulomb stress (e.g. in [18,20]), on the results is investigated in Fig. 9. and Fig. 10. Fig. 9 compares a snapshot of the results in t = 2.2. Fig. 10(a) compares the experimental and numerical, time history of maximum scour hole depth. Fig. 10 (b) compares the final (t = 3.5 s) numerical and experimental water and bed profiles. As the figures show, using Mohr–Coulomb alone results in a higher number of sediment particles to yield (to erode). This causes the model to overestimate the erosion (and consequently the scouring depth) particularly at earlier time steps when water shear stress is relatively strong. Combining the Shields yield criterion with the Mohr–Coulomb to determine whether a sediment particle moves under the applied stresses, improves the results dramatically. The mean error of the maximum scouring depth decreases from ∼0.55 to ∼0.11 by applying the additional Shields’ yield criterion. The final bed profile shows a slight underestimation of the local erosion in the proposed modification.

74

M. Khanpour et al. / Computers and Fluids 129 (2016) 67–78

Fig. 8. Continued

Fig. 9. Sample snapshot results of (a) SPH with Mohr–Coulomb yield criteria (b) SPH with Mohr–Coulomb plus Shields yield criteria, and (c) experiments, at t = 2.1 s.

3.2. Flushing test case At the initial condition, the fluid, sediment, boundary, and ghost particles were distributed uniformly with a 5 mm space (resulting in total 8872 particles). Fig. 11 shows the configuration of the SPH model setup. Initially the gate is kept closed to let the initial fluctuations damped and the hydrostatic pressure stablished in the

reservoir. The gate is then instantly removed by eliminating the ghost/boundary particles representing the gate. The peresent developed SPH rheological modeling approach which uses a combination of Mohr–Coulomb and Shields yield is applied for this test case. Fig. 12 compared sample snapshots of the numerical (SPH) and experimental results for the sediment flushing test case. Fig. 13 compares the evolution of the experimental

M. Khanpour et al. / Computers and Fluids 129 (2016) 67–78

75

Fig. 10. Comparison of experimental and SPH results using Mohr–Coulomb (MC) and Mohr–Coulomb plus Shields (MC+S) yield criteria. (a) Time history of the maximum scour hole depth (normalized using the final experimental scour depth); (b) Final water and bed profiles.

Fig. 11. SPH particle configuration for the sediment flushing test case.

and numerical water and bed surface features. Once the gate is opened, the sediment material at the frontier start to yield the critical shear stress and washed downstream creating a convex failure slope in the reservoir. The bed and free surface profiles have similar configuration in the experimental and numerical model and show a good compatibility with an insignificant error value of ∼0.04, and ∼0.16 for free surface and bed surface, respectively. Part of the discrepancy (especially at the initial stages) can be due to mechanism of the gate opening. The selection of the rheological model parameters, particularly the friction angle and Bingham viscosity, can also affect the compatibility of the results.

Fig. 14 compares the numerical (SPH) and experimental sediment flushing frontier mean slope (right before the reservoir bottom outlet). The mean slope have been calculated by matching a linear trend to the bed profile data. Fig. 14a shows, for an example time of t = 1.0 s, the numerical and experimental results have a similar slope profile, although the experimental profile shows a smoother variation. As the Fig. 14b shows that the numerical and experimental mean slopes have a very similar trend, but overall the experimental profile is slightly steeper. The results conclude that the applied rheological model using both Sheilds and MohrCoulomb yield criterion have been able to predict the overall shape and variation of flushing slope profile.

76

M. Khanpour et al. / Computers and Fluids 129 (2016) 67–78

Fig. 12. Sample snapshots of numerical and experimental results for the reservoir flushing test case.

4. Conclusion In this paper a multiphase mesh-free Lagrangian model (based on smoothed particle hydrodynamics SPH method) was developed to study the mechanism of sediment transport induced by rapid water discharge. The particular interests of this study were sediment scouring downstream of a wall-jet, and reservoir sediment flushing. The inherent ability of SPH method in dealing with the interfacial deformation and fragmentations provided the opportunity to overcome the limitation of traditional mesh-based numerical for simulation the complexities involved in the multiphase sediment-water system. The multiphase model of this study was based on the continuum description of the sediment phase, where the sediment body was treated as a viscous non-Newtonian (Bingham plastic) fluid. In addition to the commonly applied Mohr–Coulomb yielding criteria, this study employed the Shields criteria (corrected by the sediment interface slope) in combination with the boundary-layer theory to resolve the near interface incipient motion of sediment particles. Algorithms for implementation of the Bingham plastic rheological model in SPH (e.g., the method of determining the

interface sediment particles, calculating effective pressure, effective viscosity, and sediment profile slope), were also developed. To validate and evaluate the numerical results, the small-scale experiments, on the wall-jet sediment scouring and reservoir sediment flushing, were set up. The results of this study showed the ability of the model to reproduce the highly deformed water and sediment profiles. Comparing the experimental a numerical results showed a very good compatibility of sediment/water features evolutions for both of the test cases, when the proposed modification of the yield criteria is applied. The results showed the formation of a scour hole downstream of reservoir bottom outlet (for the wall-jet scouring test case) which has a similar shape and size in both numerical and experimental results. An average error of about 10.8% was observed (comparing the numerical and experimental scour depth), which is acceptable considering the complexity of the phenomena. The numerical model was also able to predict water surface wave propagation, and bed sediment deposition features with relatively good accuracy. The numerical model and experimental observation showed the formation of a convex failure slope at the frontier of the sediment flushing at the reservoir. This slope, which becomes

M. Khanpour et al. / Computers and Fluids 129 (2016) 67–78

77

Fig. 13. Evolution of experimental and numerical water and bed profile features for sediment flushing test case.

stable as goes on, has very similar features in experiments and numerical model. The simulation result of this study showed that, for test cases of this study, the features of the phenomenon can be reproduced with a good degree of accuracy using the developed SPH model. An ongoing research will work in extending the model of this study for three-dimensional cases, taking advantage of a high performance parallel algorithm. The results of this study are extendable to the sediment transport problems with similar natures. References

Fig. 14. Comparison of experimental and numerical sediment frontier (0.7≤ x ≤1.0) mean slope. (a) Bed profile and mean slope for t = 1 s; (b) Time variation of the mean slope (s).

[1] Benkhaldoun F, Sahmim S, Seaid M. Solution of the Sediment transport equations using a finite volume method based on sign matrix. SIAM J Sci Comput 2009;31(4):2866–89. [2] Lin BL, Falconer RA. Numerical modelling of three dimensional suspended sediment for estuaries and coastal waters. J Hydraul Res 1996;34(4):435–56. [3] Gharehbaghi A, Kaya B. Simulation of bed changes in rivers with finite volume method by kinematic wave model. Int J Eng Appl Sci 2011;3(3):33–46. [4] Shields A. Anwendung der Aehnlichkeitsmechanik und Turbulenz forschung auf die Geschiebewegung. Mitteilung Preussischen Versuchanstalt Wasser, Erd, Schiffbau 1936 Berlin,No.26 (in German). [5] Yalin MS. Mechanics of sediment transport. Oxford, Pergamon Press 2009; 1972. [6] Kolahdoozan M, Falconer RA. Three dimensional geo-morphological modelling of estuarine waters. Int J Sediment Res 2003;18(1):1–16. [7] Chang HK, Liou JC. Discussion of a fall velocity equation, by John P.Ahrens. J Waterway, Port, Coast, Ocean Eng 2001;127(4):250–1. [8] Cheng NS. Simplified settling velocity formula for sediment particle. J Hydraul Eng 1997;123(2):149–52. [9] Van Rijn LC. Sediment pick up function. J Hydraul Eng 1984;110(10). [10] Van Rijn LC. Sediment transport partI: bed load transport. J Hydraul Eng 1984;110(10):1431–56.

78

M. Khanpour et al. / Computers and Fluids 129 (2016) 67–78

[11] Van Rijn LC. Sediment transport partIII: bed form and alluvial roughness. J Hydraul Eng 1984;110(12):1733–54. [12] Celik I, Rodi W. Suspended sediment-transport capacity for open channel flow. J Hydraul Eng 1991;117(2):191–204. [13] Meyer-Peter E, Muller R. Formula for bed load transport. In: Proceedings of international association for hydraulic research, 2nd meeting, Stockholm; 1948. [14] Phillips BC, Sutherland AJ. Spatial lag effects in bed load sediment transport. J Hydraul Res 1989;27(1):115–33. [15] Wu W, Rodi W, Wenka T. 3D numerical modeling of flow and sediment transport in open channels. J Hydraul Eng 20 0 0;126(1):4–15. [16] Wang ZY. Experimental study on scour rate and river bed inertia. J Hydraul Res 1999;37(1):17–37. [17] Fang HW. Three dimensional calculations of flow and bed load transport in the Elberiver. Institute of hydromechanics Report no. 763. Germany: University of Karlsruhe; 20 0 0. [18] Shakibaeinia A, Jin YC. A mesh-free particle model for simulation of mobilebed dam break. Adv Water Resour 2011;34(6):794–807. [19] Shakibaeinia A, Jin YC. MPS mesh-free particle method for multiphase flows. Comput Methods Appl Mech Eng 2012;229–232:13–26. [20] Shakibaeinia A, Jin YC. Lagrangian multiphase modeling of sand discharge into still water. Adv Water Resour 2012;48:55–67. [21] Gingold RA, Monaghan JJ. Smoothed Particle Hydrodynamics: theory and application to non-spherical stars. Month Notices R Astron Soc 1977;181:375–89. [22] Lucy LB. Numerical approach to testing the fission hypothesis. Astron J 1997;82:1013–24. [23] Koshizuka S, Oka Y. Moving-particle semi-implicit method for fragmentation of incompressible fluid. Nucl Sci Eng 1996;123:421–34. [24] Morris J, Fox JP, Zhu Y. Modeling low reynolds number incompressible flows using SPH. J Comput Phys 1997;136(1):214–26. [25] Khanpour M, Zarrati AR, Kolahdoozan M. Development of a SPH model for solving basic fluid mechanics problems. Sci Ser Data Rep 2013;5(2):2–20. [26] Bachmann M, Helluy P, Jung J, Mathis H, Müller S. Random sampling remap for compressible two-phase flows. Comput Fluids 2013;86:275–83. [27] Khanpour M, Zarrati AR, Kolahdoozan M. Numerical simulation of the flow under sluice gates by SPH model. Sci Iran Trans A Civil Eng 2016;21(5):1503–14. [28] Chen W, Qiu T. Simulation of fluid flow through porous media using smoothed particle hydrodynamics method. Geo-Front 2011:4195–203. [29] Zhu Y, Fox PJ. Smoothed particle hydrodynamics model for diffusion through porous media. Transp Porous Media 2001;43(3):441–71. [30] Rooka R, Yildiza M, Dosta S. Numerical heat transfer, part b: fundamental. Int J Comput Methodol 2007;51(1):1–23.

[31] Liu MB, Feng ZM, Guo ZM. Modified SPH method for modeling explosion and impact Problems. In: APCOM & ISCM 11-14th December, Singapore; 2013. [32] Monaghan JJ, Kocharyan A. SPH simulation of multi-phase flow. J Comput Phys Commun 1995;87:225. [33] Colagrossi A, Landrini M. Numerical simulation of interfacial flows by smoothed particle hydrodynamics. J Comput Phys 2003;2:448–75. [34] Hu XY, Adams NA. A multi-phase SPH method for macroscopic and mesoscopic flows. J Comput Phys 2006;213(2):844–61. [35] Shao S. Incompressible smoothed particle hydrodynamics simulation of multifluid flows. Int J Numer Methods Fluids 2012;69(11):1715–35. [36] Xenakis AM, Lind SJ, Stansby PK, Rogers BD. An incompressible SPH scheme with improved pressure predictions for free-surface generalised Newtonian flows. J Non-Newton Fluid Mech 2015;218:1–15. [37] Manenti S, Sibilla S, Gallati M, Giordano A, Guandalini R. SPH simulation of sediment flushing induced by a rapid water flow. J Hydraul Eng 2011 Submitted January 25. [38] Shao S, Lo EYM. Incompressible SPH method for simulating Newtonian and non-Newtonian flows with a free surface. Adv Water Resour 2003;26:787– 800. [39] Fourtakas G, Rogers BD, Laurence DR. Modelling sediment resuspension in industrial tanks using SPH. La Houille Blanche 2013(2):39–45. doi:10.1051/lhb/ 2013014. [40] Ataie-Ashtiani B, Shobeyri G. Numerical simulation of landslide impulsive waves by incompressible smoothed particle hydrodynamics. Int J Numer Methods Fluids 2008;56:209–32. [41] Zanganeh M, Yeganeh-Bakhtiary A, Abd Wahab KA. Lagrangian coupling twophase flow model to simulate current-induced scour beneath marine pipelines. Appl Ocean Res 2012;38:64–73. [42] Monaghan JJ. Simulating free surface flow with SPH. J Comput Phys 1994;110:399–406. [43] Shakibaeinia A, Jin Y. MPS-based mesh-free particle method for modeling open-channel flows. J Hydraul Eng 2011;137(11):1375–84. [44] Liu MB, Liu GR. Smoothed particle hydrodynamics (sph): an overview and recent developments. Arch Comput Methods Eng 2010;17:25–76. [45] Shahriari S, Hassan IG, Kadem L. Modeling unsteady flow characteristics using smoothed particle hydrodynamics. Appl Math Model 2013;37(3):1431–50. [46] Shirkhani H, Khanpour M, Zarrati AR. Development of mesh-free particle technique for simulation of non-hydrostatic free surface flow. Comput Fluids 2014;105:8–15. [47] Van Rijn LC. Principles of sediment transport in rivers, estuaries, and coastal seas. The Netherlands: Aqua Publications, Blokzijl; 1993.