Simulations of Local Scour Around a Cylindrical Bridge Pier And a Semicircular Abutment Using Unsteady k-ε Model Combined with σ-Grid

Shi Liu1, *, Yi Yang1, Xiaobo Wu2
1 Electric Power Research Institute of Guangdong Power Grid Co., Ltd., Guangzhou, Guangdong, China
2 Wuhan DaCheng Construction Consulting Co., Ltd., Wuhan, Hubei, China.

Article Metrics

CrossRef Citations:
Total Statistics:

Full-Text HTML Views: 1269
Abstract HTML Views: 508
PDF Downloads: 442
ePub Downloads: 227
Total Views/Downloads: 2446
Unique Statistics:

Full-Text HTML Views: 603
Abstract HTML Views: 272
PDF Downloads: 284
ePub Downloads: 143
Total Views/Downloads: 1302

© 2017 Liu et al.

open-access license: This is an open access article distributed under the terms of the Creative Commons Attribution 4.0 International Public License (CC-BY 4.0), a copy of which is available at: This license permits unrestricted use, distribution, and reproduction in any medium, provided the original author and source are credited.

* Address to correspondence to this author at Electric Power Research Institute of Guangdong Power Grid Co., Ltd., Guangzhou, Guangdong, China; Tel: (+86)18627931026; E-mail:



Bridge scour is one of the major causes of bridge failure. Prediction of the maximum depth and shape of local scour plays an important role in bridge design and maintenance. In this paper, k-ε turbulent model combined with wall function was employed and complex flow fields are exhibited. The bed load transport model summarized by Qian and Wan [1] was applied to evaluate the development of local scour. Besides, σ-grid module was developed and embedded into the commercial solver FLUENT to fit the change of channel-bed.


In this module, the change of the elevation of the channel bed is calculated through using user defined functions(UDF), and the motion of the grid system is realized based on a program developed through C++ language, which extends the applications of FLUENT. The evolution of local scour hole for two cases, i.e., local scour around a cylindrical pier and a semicircular abutment, was simulated.

Results and Conclusion:

The depth and the shape of local scour as well as the flow fields were predicted. Numerical result conforms well to the experimental measurement. Especially, it provides fairly reasonable prediction on the key issue of the maximum scour depth. The satisfactory agreement validates the numerical method developed in the present study. In addition to the validation purpose, the different performance of this method for predicting local scour around the cylindrical pier and the semicircular abutment was discussed.

Keywords: local scour, k-ε model, bridge pier, abutment, σ-grid, sediment transportation, FLUENT.


Sediments on the channel bed interact with energetic coherent vortex structures leading to local scour which could undermine the structural reliability of embedded hydraulic structures and endanger the channel bed stability. For bridge engineering, local scour around bridge piers or abutments is commonly recognized as one of the major causes of bridge failures at river crossing. Local scour has led to massive lose on life and economy, and resulted in serious impact on local transportation as pointed out by Shepherd [2]. Therefore, the prediction of local scour especially the maximum depth of scour hole, is critical for bridge design, maintenance and evaluation. Over the past decades, many experiments have been conducted to study this issue as Yanmaz [3], Ettema [4], and Sheppard [5, 6]. However, the detailed vortex structure is difficult to be revealed based on experimental measurements. In addition, there are also several scour equations based on laboratory data but they did not always produce reasonable results according to Chavan [7]. In this context, Gaudio et al. [8, 9] showed that due to different mathematical structure of empirical formulas, different results may be obtained through using scour depth predictive equations. Numerical simulation is another approach to examine the flow structure in the hole of local scour or the development of local scour. In the present study, the numerical approach is adopted. The numerical simulation concerning local scour could be grouped in two types. One type is the fixed bed simulation, the other is the movable bed simulation. In the following, we briefly review these two types of simulations.

Due to the fact that there is no variation of the elevation of river bed, the computation for flow fields over the fixed bed is not as useful as the live-bed case. Therefore, we can use more accurate turbulent models to study this issue, such as large eddy simulation (LES). Zhao and Aode [10] applied LES turbulent model to examine the flow fields and the scouring mechanism around a circular pier. They found that at the downstream of the pier, when the scouring depth increases, the turbulent intensity as well as the pressure fluctuation increases. Besides, the bed shear stress at the upstream of the pier is primarily influenced by the horseshoe vortices. Pasiok and Stilger [11] investigated the mass particle trajectories around a pier through large eddy simulation and concluded that the particle trajectories depend on intensity of vertical vortices which is the main driving force moving the particles out of the scour hole. Bressan et al. [12] applied large eddy simulation to examine the turbulent flow around a scoured bridge abutment and concluded that the dynamics of the primary vertical structure can induce large shear stress on the bed, leading to sediment transport to be highly intermittent. Kirkil and Constantinescu [13, 14] studied the flow and turbulence structure around an in-stream rectangular cylinder and a spur dike with scour hole on the basis of detached eddy simulation (DES). In the studies, they highlighted the need to understand and differentiate the mechanisms causing scour around bluff body obstacles of different shapes and at different stages of the scour process.

For the movable bed simulations, currently, there are three approaches. The first one is to consider the sediment and water as two phases. The second approach is based on the bed shear stress and the sediment transportation model and the last is based on smoothed particle hydrodynamics (SPH) model.

For the first approach, Richardson et al. [15] applied the two-phase flow modeling techniques and simulated the flow structures around a bridge pier with and without the scour hole. Compared the simulated results with the experiment ones, they found that the FLOW3D hydrodynamic model well simulates the complex flow patterns around the bridge pier. Huang et al. [16] applied two-phase Eulerian model to simulate the local scour progress. They examined the scale effects on sediment scour and found that ignoring Reynolds similarity in physical modeling may result in errors for scour in large bridge piers. However, two-phase model will increase the number of governing equations, thus definitely increasing the computing time.

For the second approach, Zhao et al. [17, 18] employed a three-dimensional finite element numerical model to simulate the local scour around submerged vertical cylinders. k-ω turbulence closure was applied to solve the URANS equations and found that the horseshoe vortex combined with the other vortices governs the scouring process. Abdallah et al. [19] applied SSIIM program and solved Reynolds-stress term by k-ε model, indicating that the curvature shape of bridge abutment could reduce the local scour depth by more than 95%. Ehteram and Meymand [20] selected SSIM2 code to model the scour depth at side piers of the bridge. They found that the SSIM2 model slightly underestimates the scour depth but overall the results from SSIM2 are acceptable. Basser et al. [21] applied SSIIM 2.0 to model the scour around a rectangular abutment. SSIIM 2.0 employs k-ε turbulence model with some RNG extensions for predicting the turbulence. It could be concluded that SSIIM 2.0 numerical modeling is capable of simulating scour around the rectangular abutment. Xiong et al. [22] used commercial software FLUENT together with the dynamic mesh available in FLUENT to simulate the local scour around bridge piers. However, the tetrahedral mesh with quite large size was applied near the bridge piers and the channel bed which will influence the accuracy of the prediction of local scour depth.

Mirmohammad and Ketabdari [23], Ran et al. [24] presented the method simulating local scour by SPH model. However, up to now, in the most commercial solvers, the SPH model has not been integrated, limiting the applications.

In the present study, the commercial code FLUENT (Fluent 6.3, Fluent Inc.) was extended by user defined functions (UDFS) to make the function of modeling local scour with structured mesh available in it. The block structured mesh will be used with very fine mesh close to the pier and the channel bed. The second numerical approach such as solving the sediment continuity equation, was employed and the change of topography was fitted by σ grid technique. Two cases (cylindrical pier and semicircular abutment) will be examined to verify the model developed. In addition, the different scour patterns between them are also discussed.


The governing equations are introduced in section 2.1. Section 2.2 and section 2.3 present the sediment transportation model and the σ grid respectively, followed by the presentation of the boundary conditions used in this study.

2.1. Governing Equations

The turbulent flow is simulated by solving the Unsteady Reynolds-Averaged Navier-Stokes (RANS) equations with a k-ε turbulence model. The Arbitrary Lagrangian Eulerian (ALE) scheme was employed for solving the RANS in a computational domain with continuously developing bed boundary. In the ALE scheme, the position of each mesh node moves and convection terms are modified to consider the effects of the mesh moving speed. According to the FLUENT 6.3 User’s Guide [25] the URANS are expressed as:


in which, xi (i=1, 2, and 3) is the coordinate along axis in Cartesian coordinate system, ui is the instantaneous velocity component in xi direction; ρ is the density of clear water; t is time; ujp is the velocity of the movement of the computational mesh; ν is the kinematic viscosity; τij is the Reynolds stress defined as τij = νt(∂ui/∂xj + ∂uj/∂xi)−(2/3)ij; νt is the turbulent viscosity; k is the turbulent kinematic energy and defined as k = 1/2u'iu'i; u'i is velocity fluctuation in xi direction. The standard k − ε turbulent model is employed to close the RANS. For the standard k − ε turbulent model, turbulent kinetic energy k and turbulent dissipation ε rate are obtained based on the following transport equations:


where represents the generation of turbulence kinetic energy due to the mean velocity gradients; μt = ρCμk2 / ε; C1ε, C2ε, Cμ, σk, σε are the empirical constants, with value of C1ε = 1.44, C2ε = 1.92, Cμ = 0.09, σk = 1.0, σε = 1.3 in this study.

In standard k − ε turbulent model, the flow is assumed to be fully developed turbulent flow. In order to avoid fine grids among sub-layer and save computational time, wall function is applied to obtain wall sheer stress. It can be done by using:


In which subscript “p” stands for the center of control volume adjacent to wall; Up is the velocity component parallel to the wall; δp is the distance from the center of control volume to the wall; Γwall is the dissipation constant and can be obtained from the following equation:


where is turbulence kinetic energy at point “P”; κ is the Karman constant. Gaudio etc. [26] argued that the Kármán constant is however nonuniversal in flows over mobile sediment beds. Besides, Foken[27], Hogstrom[28, 29] report values of κ between 0.35 and 0.42. Following the study by Fuhrman[30]. 0.4 is chosen for the Karman constant. E = eκ(B−ΔB) is the roughness empirical constant; B= 5.2 is additive constant; ΔB is the roughness function which reads as follows:


In which ks+ = u*ks / ν is the effective roughness height of the wall boundary or particle Reynolds number playing an important role in quantifying the influence of roughness on the bed and solid boundaries. For the stationary flat bed with uniform grain size, ks should be equal to the median grain size, d50. At present, a value of ks = d50 for the bed is employed.

2.2. Sediment Transportation Model

The scour problem is due to the movement of the channel bed sediments. Bed materials that are transported along the bed are called bed load. Bed load moves by rolling, sliding, and hopping. However, if particles settle down slowly enough and are carried in current either never touching the bed or only intermittently touching it, they are treated as suspended load. The scouring type, i.e., bed load and suspended load, depends on the size of particles and flow velocity. If the average diameter of the particles is larger than the critical diameter, sediments will be transported in the form of bed load. Otherwise, suspended load will be the main transportation form. The critical diameter of particles can be obtained by the equation provided by Qian and Wan [1] as:


where U is the mean approaching flow velocity; d is the diameter of sediments. In this paper, the mean approaching flow velocities for the bridge pier case and abutment case are 0.25m s-1 and 0.294m s-1, with sediments average diameter of 0.385mm and 0.52mm respectively, which is significantly larger than the corresponding values of U2 / (360⋅g). As a result, bed load is considered as the main transportation form. Such physical processes are governed by the sediment continuity equation combined with the bed-load transport equation. We applied herein the bed morphodynamic model summarized by Qian and Wan [1]. According to sediment transport balance, the variation of bed elevation can be derived as:


where γ is the sediments porosity; h is the elevation of the channel bed; and qbi represents the bed load flux vector in the xi direction. If we know the bed load flux at each time step, the variation of the channel bed elevation could be calculated using Eq.(9). Explicit Euler scheme is chosen to solve Eq(9), indicating that the renewed elevation of channel bed is calculated using the information of the flow fields at the previous time step with the consideration of the easy of programming. If implicit scheme is applied, during iterations not only the flow fields will be renewed by the grid coordinates which bringing much more difficulties. In the future, we will attempt to use more accurate implicit schemes to solve the N-S equation as well as the sediment transportation equation at the same time. In fact, an accurate prediction of the bed load flux plays a critical role in correct estimation of scour depth variation. Recently, Dodaro et al. [31, 32] modified Einstein bed load equation to predict local scour downstream of a rigid bed. In the present study, bed load flux was estimated as follows. Considering the slope of channel bed, Eq.(10) is adopted to determine qbi following the study by Qian and Wu[1]:


in which τi is the bed shear stress in xi direction; |τ| is the magnitude of the bed shear stress; C = 1.5 an empirical constant following the study by Qian and Wu[1]; h is the bed elevation; q'b is determined from:


in which diis the diameter sediment whose weight is less than i% of all sediment; ξ = (ρsρ) / ρ is dimensionless density; ρs is the sand density; g is the gravitational acceleration 9.8m s-2; D* = d50 [(s − 1)g / ν2]1/3 is the dimensionless particle parameter; s = ρs/ ρ is the relative density; T is the non-dimensional excess bed shear stress; and Tm is the modified non-dimensional excess bed shear stress; T and Tm are expressed as:


where λ = e(0.45α + 0.2β) is correction factor; α = 1 − σg; β = d50 / d90d10 / d50; σg = (d84 / d50 + d50 / d16)/ 2; τ is the shear stress u* is the shear velocity; τb,cr is the critical bed-shear stress at which sediment starts to move. For the local scour problem, Qian and Wan [1] presented a formulation of critical sheer stress considering the effects from the bed slope:


where φ is the angle of repose for sediments; CL donates the lift coefficients for the particle; CD is the drag coefficients for the particles; CL/CD=0.85, as indicated by Ye and Wen [33]; α is the angle between the flow direction and the slope direction projected onto a horizontal plane; θ is the angle between the slope and horizontal plane, as shown in Fig. (1). τ'b,cr is critical bed-shear stress for sediments on horizontal bed and can be obtained by:

Fig. (1). Forces acting on a sediment particle on an inclined plane.


in which θcr is the Shields number, which can be obtained from the shields curve. The sediment transportation model presented above is coded using user defined functions in FLUENT.

2.3. σ-Grid

For the live-bed calculations, the topography of the channel bed changes with time. Therefore, in order to fit channel bed at each time step a novel adaptive grid method, σ adaptive grid is applied. The basic philosophy of the method can be summarized as follows. In the σ adaptive grid system, the ratios of the neighboring grid sizes hold as constants. An example of two dimensional σ adaptive grid system can be found in Fig. (2a and 2b). Fig. (2a) shows the grid at the stage without the variation of the elevation of the channel bed and Fig. (2b) is with the variation. It could be found that the total grid number and the horizontal grid size do not change before and after scouring. What changed is only the grid coordinates in the vertical direction but the ratios of the neighboring grids hold as constants. In this way the variation of the channel bed topography is fitted quite well. The σ adaptive grid is not available in FLUENT. Therefore, with the intent to apply σ adaptive grid, we extended FLUENT using language C++ which could control the running of FLUENT and read or change the node coordinates at each time step. The time step size is set as 0.001s. Several time step sizes have been examined and found that the numerical results will be independent with the time step size when it is smaller than 0.001s.

Fig. (2). Two-dimensional example of the appliction of σ grid fitting the variation of the channel bed, (a) before scouring and (b) after scouring.

The procedure in the calculations was as follows:

Step 1. generate the mesh using Gambit;

Step 2. calculate the flow;

Step 3. calculate the sediment transport;

Step 4. calculate the variation of channel bed elevation;

Step 5. update the channel bed using σ-grid;

Step 6. return to Step 2.

2.4. Boundary Conditions

The computational domain includes boundary conditions of inflow, outflow, symmetry boundaries and walls.

2.4.1. Inflow and Outflow

For each inlet boundary, Dirichlet boundary condition as has been introduced in detail by Ferziger (2002) is given and the transverse as well as the vertical velocity components are set to be zero. In each case, a separate calculation of flow in a straight channel of a length of 20D (diameter of the pier or abutment) is carried out whose inlet velocity is determined by U = Q / A, where Q is the water discharge and A is the wet cross-section area, and turbulent kinetic energy k as well as turbulent dissipation rate ε at the inlet are 0 respectively. The velocities as well as the turbulent quantities at the outlet boundary of the straight channel are used as the inlet boundary for scour simulation. This is a reasonable way to produce a stable boundary layer flow condition constituent to the rough riverbed and the turbulent model used. At the outlet boundary, the Neumann conditions were applied. The gradients of the velocity and the turbulent quantities ( / n) in the flow directions are set to zero and the reference pressure of zero is provided. We did not specify the velocity values at the outlet. The flow depth is constant for all of computational domain. Multiphase technique could be used to calculate the free surface. However, due to the limited computational resource, we finally decide to use the symmetry boundary condition for the modelling of free surface, so the water depth is a constant.

2.4.2. Symmetry and Wall Boundaries

At two side boundaries for the pier case, one side boundary for the abutment case as well as the top boundaries for both two cases (water surface boundary), as shown in Fig. (3), the velocity component and gradients of other dependent variables in normal direction are set as zero. At the bed for the pier case and the abutment cases as well as the side attaching the abutment, the standard wall function boundary condition is adopted. The logarithmic law considering the roughness on the wall was used and depicted as Eq.(7).

Fig. (3). Computational domain and boundary conditions for the numerical simulation of the scour around (a) a cylindrical pier and (b) a semicircular abutment.


3.1. Numerical Models

Two experimental models, i.e. cylindrical bridge pier by Melville [35] and the semicircular abutment by Dey and Karim[36], were numerically modeled.

3.1.1. Cylindrical Bridge Pier

To validate the effectiveness of the present method, numerical simulations of flow and bed deformation around a cylindrical pier are carried out firstly. The simulation condition is almost the same as that performed in laboratory experiments reported by Melville [35]. The tests were conducted in a water flume of 19m long and 45.6cm wide. A cylinder with diameter, D, of 5.08cm was simulated as a bridge pier model. The bed material involved was relatively uniform sized sand with a median grain size of 0.385mm. The approaching mean flow velocity is 0.25m s-1. Water is discharged at the rate of 0.01712m3 s-1 with water depth of 0.15m. Fig. (4a) shows the sketch of the experiment setup. For the computational domain in the numerical simulation, inlet locates at a distance of 7D upstream of the cylinder whose diameter is the same as that in the experiment by Melville (1975). The outlet locates at a distance of 10D downstream the cylinder. The computational domain is 0.15m thick and 45.6cm wide which are identical with those in experiment. The origin point of the domain locates at the center of the cylinder attaching the original channel bed.

For the abtment case, the experiment by Dey and Karim [36] is modeled. That experiment was conducted in a flume with length of 20m, width of 0.9m and depth of 0.2m. A semicircular abutment model with a diameter of 0.26m is located in this flume in which the sand with uniform diameter of 0.52mm is filled to model the channel bed, as shown in Fig. (4b). The water depth at the inlet is 0.2 and the mean flow speed is 0.294 m s-1. As a result, the flow rate is 0.0524 m3 s-1. For the computational domain in the numerical simulation, inlet locates at a distance of 7D upstream of the semicircular abutment. The outlet locates at a distance of 10D downstream. The computational domain is 0.2m thick and 0.9m wide which are identical with those in experiment, as shown in Fig. (4b). And the origin point of the domain locates at the center of the semicircular abutment attaching the original channel bed.

Fig. (4). Schematic plot of test arrangement of (a) the cylindrical pier and (b) the semicircular abutment.

3.2. Mesh System

The generation of an appropriate grid has been proved to have a major effect upon model prediction and ideal grid would produce a flow field which is independent of grid dimensions, with as low grid density as possible to save computational time. The computational domain is discretized by multi-blocked structure hexahedral grids. In the region close the walls of the structures (pier and abutment) and the channel bed, the flow is very complex. Therefore, the grids are refined. Near the walls and the ground the first grid layer size is set as 0.2mm(fine), 0.5mm(medium), and 1mm(coarse) to check if the numerical simulation is grid independent through comparing the flow fields of the fixed bed situations obtained from each grid system. It is found that the fine and medium grid systems give almost the same results, indicating that the grid independence requirement is reached at the medium grid stage. In the present study, the results of the fine grid system will be presented. The grid growing ratio is 1.1 to avoid the numerical errors resulting from the sudden change of the mesh. The largest grid size is limited as 2cm. Totally,105948 grids are generated for the calculation of the pier and 201526 for the abutment. The detailed picture of the grid distribution could be found in Fig. (5).

Fig. (5). Detailed of the mesh arrangement for (a) the cylindrical pier and (b) the semicircular abutment.

3.3. Solving Schemes

The governing equations are discretized on collocated grid system through using the finite volume method. The second order upwind scheme is used for the convective and viscosity term. The central difference scheme is used for the diffusive term. SIMPLEC algorithm is employed as introduced by Ferziger and Peric [34]. Time step size is 0.001s. The computation was stopped after 30min for the pier case and 60min for the abutment case because the maximum change of the scour hole in between 1min is less than 0.01mm. The computational time is huge. For single simulation, it will cost 7 days on the cluster with 8 nodes and each node consists of 8cores. The floating point speed of each core is 2.4GHz.


4.1. Cylindrical Bridge Pier

Fig. (6a, b and c) show the grid systems of the bridge pier with scour at 10min, 20min and 30min respectively. In Fig. (6), we would like to show the change of channel bed more clearly; therefore, only 1/4 of total grids are shown. It can be found that σ-grid fits the complicate terrain successfully. At 30min, the topography is most complex and the gradient of the elevation of channel bed gives the maximum value, however, the grid quality is still very good with the maximum skew ratio of only 0.1.

Fig. (6). Grids in around the cylindrical pier at (a) 10 min, (b) 20 min, and (c) 30 min.

Fig. (7) shows the contour of the scour elevation and the corresponding experimental results. The difference between the computed maximum equilibrium scour depth (approximately 5cm) and the corresponding measured values of 4.2cm is within 25% . A significant qualitative difference does exist between the experiment and the simulation. The former shows a maximum scour region wrapping around the front of the pier and extending up until about 80o angle from the streamwise central line. Nevertheless, the simulation shows that the maximum scour is confined in a small region at the side of the cylinder. This discrepancy should be due to the inherent inability of URANS models to resolve the dynamics of the horseshoe vortex system at the nose of the cylinder which has historically been considered as the main driving force of scour in this region. The horseshoe vortex system will be discussed later in more detail to demonstrate the mechanism of the difference between the simulation and experiment.

Fig. (7). Local scour shape at t=60min. (a) numerical results and (b) experimental results. (cm).

Fig. (8a and b) show comparisons of the distribution of the flow speed on the slice of y=0 between the experiment measurements and the numerical simulations at the equilibrium stage. The data from simulation is time averaged. After the scour reached the equilibrium stage, we stopped the topography calculation but continued the simulation for another 2min. The flow fields at each time step during this 2min are sampled to obtain the time averaged flow fields at the equilibrium stage. The speed is normalized the same as the experiment by the mean approach flow velocity that was 0.25m s-1. In both the simulation and the experiment, it can be found that a strong downward flow developed along the pier face produced rather large velocity components near the bed, which distorted the profiles in the vicinity of the cylindrical pier. Flow separates at the nose of the scour hole and reattaches at the front of the pier, forming a horseshoe vortex, which can be clearly identified from both the simulation and the experiment.

The horseshoe vortex induces high bed shear stress, leading to the suspension of the sediment in front of the cylinder. The suspended sediment is washed away to the downstream of the cylinder by this vortex. Nevertheless, in comparison with the experiment, we can find that the horseshoe vortex does not fully develop and the dimension of the vortex is weakened. In the experiment, the horseshoe vortex reaches to the location just beneath the cylinder, and the vertical size of the vortex is about twice as large as that in the simulation, implying that the URANS model could not resolve the dynamics of the energetic horseshoe system well at the upstream junction of the cylinder with the bed. More sophisticate turbulent model, such as LES (large eddy simulation), and sediment transportation model should be applied in the future.

In fact, LES model has a great ability to resolve vortices. Recently, Calomino et al. [37] successfully simulated small vortices within the macro roughness elements in a corrugated pipe by means of LES model. As a result, it is believed that LES has an acceptable capability to resolve horseshoe vortex around the pier. However, LES model requires larger grid density and smaller time step size which indicates the computational time will increase significantly. We have pointed out that using the 8nodes×8cores cluster combined with the URANS model the computation will last for 7 days. Thus, if the LES model is used, based on the computational resources applicable, now the computation time will increase to at least several months. Therefore, we can conclude that for the numerical prediction of the scour, the URANS model is currently more suitable than LES from engineering application point of view.

Fig. (8). Velocity contours at t=60min in front of the cylindrical pier. (a) numerical results and (b) experimental results. (m s-1).

4.2. Semicircular Bridge Abutment

Admittedly, it is the first time to numerically simulate the local scour around a semicircular bridge abutment. Three dimensional view of the simulated scour hole around a semicircular abutment is shown in Fig. (9) with positive x being the inflow direction. It can be found that the scour depth increases as closing the abutment and shows the maximum depth at the lee side, after which the scour depth decreases, proving a reasonable shape of local scour.

Fig. (9). 3-D view of the local scour around the semicircular abutment from numerical simulation at t=60min.

Six slices with azimuthal angle of 10°(A-A), 30°(B-B), 45°(C-C), 90°(D-D), 120°(E-E) and 170°(F-F) same as those in the experiment by Dey and Barbhuiya(2005) are shown in Fig. (10). The profiles of the scour and the flow fields on these six slices will be discussed. In the following discussion, r is the radiation distance from center of the abutment, θ is the azimuthal angle. Same as experiment the vertical coordinate and the radial distance will be normalized by l=D/2.

Fig. (10). Six vertical slices for plotting local scour profile and the flow fields around the semicircular abutment.

Different from the previous cylinder pier case, horseshoe vortex system around the semicircular abutment is not developed very well as has been pointed out by Gogus [38]. Therefore, it is reasonable to anticipate that scour in this case should primarily be driven by the bed shear-layers that sheds from the nose of the abutment and locally increases at the lee side.

Besides, it is also anticipatable that the URANS model should be more appropriate for the prediction of the local scour around the abutment. Fig. (11) shows the comparison of the scour hole at the equilibrium stage, from which we can see the simulated results give quite good agreement with experiment. In the simulation the maximum depth of the scour hole, dm, reaches to 153mm which is only 9mm lower than the experiment. The bed elevation patterns are very comparable based on the measurements. In fact, for this case, the simulation captures correctly the location of maximum scour patterns and the overall shapes of the bed elevation contours.

Fig. (11). Top view of local scour around the semicircular abutment from (a) Experimental results and (b) Numerical results, at t=60min. (mm).

In addition, flows on two horizontal slices with z = -0.5dm and z=0 (dm: the maximum depth of scour) were also measured. The comparison about the flow on these horizontal slices is also carried out as shown in Fig. (12) and Fig. (13). At z = -0.5dm the scour hole in experiment is measured which starts from 2.4l and ends at -2.4l which is close to those in the simulation whose scour hole starts and ends at 2.2l and -2.2l respectively. In the spanwise direction, the scour size reaches to around 2.0 in both of the simulation and the experiment. At the height of z=0, the size of the scour hole is underestimated based on numerical simulation in the streamwise direction but agrees well in the spanwise direction.

Fig. (12). Velocity vectors at t=60min on the horizontal cross section of z = -0.5dm. (a) experiment and (b) numerical simulation.

Fig. (13). Velocity vectors at t=60min on the horizontal cross section of z=0. (a) experiment and (b) numerical simulation.

Fig. (14). Velocity vectors at t=60min on vertical cross sections, (a) θ=10o in experiment, (a’) θ=10o in simulation, (b) θ=30o in experiment, (b’) θ=30o in simulation, (c) θ=45o in experiment, (c’) θ=45o in simulation, (d) θ=90o in experiment, (d’) θ=90o in simulation, (f) θ=170o in experiment, (f’) θ=170o in simulation.

The profiles of local scour on the six vertical slices are also plotted with the vectors and contour of fluid speed superimposed as shown in Fig. (14) and Fig. (15) respectively. It is clear that the shape of the scour hole at each vertical slice in simulation is comparable with that in experiment. In general, the distributions of the speed give satisfactory agreement between simulation and experiment. There are some differences been the simulation and experiment, which may be due to the fact that when the scour develops, the mesh quality will become lower, especially near the bed, because the grid will be not normal to the bed and the angle decrease to a small value at the locations in which the terrain shows large gradient. This should be the most important reason of the errors in the numerical simulation, further validating the numerical method proposed in the present study.

Fig. (15). Velocity contours at t=60min on vertical cross sections, (a) θ=10o in experiment, (a’) θ=10o in simulation, (b) θ=30o in experiment, (b’) θ=30o in simulation, (c) θ=45o in experiment, (c’) θ=45o in simulation, (d) θ=90o in experiment, (d’) θ=90o in simulation, (f) θ=170o in experiment, (f’) θ=170o in simulation.

Fig. (16a and b) are the shear stress distribution on the channel bed before and after the formation of local scour. We can see that at the equilibrium stage the shear stress around the abutment decreases significantly and are lower than the critical shear stresses.

Fig. (16). Bed shear stress (a) before scouring and (b) at t=60min from numerical simulation. (Pa).


To conclude, in this study, a method combining the k-ε turbulence model with the sediment transportation model is applied and σ adaptive grid module is developed in commercial solver FLUENT to fit the variation of the channel bed. The formation of local scour is numerically simulated. This study aims at verifying the proposed numerical method. The shape of local scour at t=60min and the details of the flow fields after the formation of local scour is examined and compared with the experiment. We find that the proposed method successfully reproduces the development of local scour and the maximum depth of the scour hole which is the critical parameter in the design of the bridge is well predicted, indicating that the k-ε model together with the sediment transportation model is applicable to estimate the local scour for the bridge engineering. Apart from the validation purpose, we also found that the prediction for the abutment is much better than that for the pier. This is considered to be due to the fact that the strong horseshoe vortex is the main driven force for the scour in front of the pier but could not be well reproduced by URANS model. Further researches using LES turbulence model combined with more advanced sediment transportation model should be conducted.


The following symbols are used in this paper:

CL  = lift coefficients for the particle;
CD  = drag coefficients for the particles;
d  = diameter of sediments. In this paper;
di  = diameter sediment whose weight is larger than i% of all sediment;
D*  = dimensionless particle parameter;
E  = roughness empirical constant;
g  = gravitational acceleration;
h  = elevation of the channel bed;
k  = turbulent kinematic energy;
kp  = turbulence kinetic energy at point “P”;
qbi  = bed load flux vector in the xi direction;
s  = relative density;
t  = time;
T  = non-dimensional excess bed shear stress;
Tm  = modified non-dimensional excess bed shear stress;
u*  = shear velocity;
ui  = instantaneous velocity component in xi direction;
ujp  = velocity of the movement of the computational mesh;
u'i  = velocity fluctuation in xi direction;
U  = mean approaching flow velocity;
Up  = velocity component parallel to the wall;
xi  = the coordinate along axis in Cartesian coordinate system;
α  = angle between the flow direction and the slope direction projected onto a horizontal plane;
γ  = sediments porosity;
Γwall  = dissipation constant;
δp  = distance from the center of control volume to the wall;
θ  = angle between the slope and horizontal plane;
θcr  = Shields number;
κ  = Karman constant;
ν  = kinematic viscosity;
νt  = turbulent viscosity;
ξ  = dimensionless density;
ρ  = density of clear water;
ρs  = sand density;
τ  = shear stress
|τ|  = magnitude of the bed shear stress;
τb,cr  = critical bed-shear stress at which sediment starts to move;
τ'b,cr  = critical bed-shear stress for sediments on horizontal bed;
τi  = bed shear stress in xi direction;
τij  = Reynolds stress;
φ  = angle of repose for sediments.


Not applicable.


Not applicable.


The authors declare no conflict of interest, financial or otherwise.


Declared none.


[1] Qian N, Wan Z, Armenio V. Sediment transport mechanics. Science Press, Beijing, 2003 (in Chinese)
[2] Shepherd R, Frost JD. Failures in Civil Engineering: Structural, Foundation and Geoenvironmental Case Studies. Am Soc Civ Eng 1995.
[3] Yanmaz AM, Altinbilek HD. Study of time-dependent local scour around bridge piers. J Hydraul Eng 1991; 117(10): 1247-68.
[4] Ettema R, Melville BW, Barkdoll B. Scale effect in pier-scour experiments. J Hydraul Eng 1998; 125(8): 894-5.
[5] Sheppard DM, Glasser T, Odeh M. Large scale clear-water local pier scour experiments. J Hydraul Eng 2004; 130(10): 957-63.
[6] Sheppard DM, Miller W. Live-Bed Local Pier Scour Experiments. J Hydraul Eng 2006; 132(7): 635-42.
[7] Chavan R, Mukhopadhyay N, Singh A, Kumar B. Comparison of existing equation for scour at bridge piers using laboratory and field data. International Conference on Hydraulics, Water Resour Coastal Eng.
[8] Gaudio R, Grimaldi C, Tafarojnoruz A, Calomino F. Comparison of formulae for the prediction of scour depth at piers. Proc 1st European IAHR congress. Edinburgh, UK. 2010.
[9] Gaudio R, Tafarojnoruz A, De Bartolo S. Sensitivity analysis of bridge pier scour depth predictive formulae. J Hydroinform 2013; 15(3): 939-51.
[10] Zhao W, Huhe A. Large-eddy simulation of three-dimensional turbulent flow around a circular pier. J Hydrodyn Ser: b 2006; 18(6): 765-72.
[11] Pasiok R, Stilger-Szydło E. Sediment particles and turbulent flow simulation around bridge piers. Arch Civ Mech Eng 2010; 10(2): 67-79.
[12] Bressan F, Ballio F, Armenio V. LES of turbulence around a scoured bridge abutment. ERCOFTAC Series. vol.15, pp. 251–256, 2011.
[13] Kirkil G, Constantinescu G. Flow and turbulence structure around an in-stream rectangular cylinder with scour hole. Water Resour Res 2010; 46(11): 2387-92.
[14] Kirkil G, Constantinescu G. Flow and turbulence structure around a spur dike in a channel with a large scour hole. Water Resour Res 2011; 47(12): 1-19.
[15] Richardson JE, Panchang VG. Three-dimensional simulation of scour-inducing flow at bridge piers. J Hydraul Eng 1998; 124(5): 530-40.
[16] Huang W, Yang Q, Xiao H. CFD modeling of scale effects on turbulence flow and scour around bridge piers. Comput Fluids 2009; 38: 1050-8.
[17] Zhao M, Cheng L, Zang Z. Experimental and numerical investigation of local scour around a submerged vertical circular cylinder in steady currents. Coast Eng 2010; 57(8): 709-21.
[18] Zhao M, Vaidya S, Zhang Q, Cheng L. Local scour around two pipelines in tandem in steady current. Coast Eng 2015; 98: 1-15.
[19] Mohamed A, Nasr-Allah H, Abdel-Aal M, Awad S. Investigating the effect of curved shape of bridge abutment provided with collar on local scour, experimentally and numerically. Ain Shams Eng J 2015; 6: 403-11.
[20] Ehteram M, Meymand AM. Numerical modeling of scour depth at side piers of the bridge. J Comput Appl Math 2015; 280: 68-79.
[21] Basser H, Cheraghi R, Karami H, et al. Modeling sediment transport around a rectangular bridge abutment. Environ Fluid Mech 2015; 15: 1105-14.
[22] Xiong W, Cai C, Kong B, Kong X. CFD Simulations and analyses for bridge-Scour development using a dynamic-mesh updating technique. J Comput Civ Eng 2014; 30
[23] Mirmohammadi A, Ketabdari M. Numerical simulation of wave scouring beneath marine pipeline using smoothed particle hydrodynamics. Int J Sediment Res 2011; 26(3): 331-42.
[24] Ran Q, Tong J, Shao S, Fu X, Xu Y. Incompressible SPH scour model for movable bed dam break flows. Adv Water Resour 2015; 82: 39-50.
[25] FLUENT . User’s Guide 63. USA: Fluent Inc. 2006.
[26] Gaudio R, Miglio R, Dey S. Nonuniversality of von Kármán’s κ in fluvial streams. J, Hydraul Res, International Association for Hydraulic Research 2010; 48(5): 658-63. [IAHR].
[27] Foken T. 50 years of the Monin-Obukhov similarity theory. Boundary-Layer Meteorol 2006; 119: 431-47.
[28] Hogstrom U. Non-dimensional wind and temperature profiles in the atmospheric surface layer-a re-evaluation. Boundary-Layer Meteorol 1988; 42: 55-78.
[29] Hogstrom U. Review of some basic characteristics of the atmospheric surface layer. Boundary-Layer Meteorol 1996; 78: 215-46.
[30] Fuhrman DR, Baykal C, Sumer BM, Jacobsen NG, Fredsøe J. Numerical simulation of wave-induced scour and backfilling processes beneath submarine pipelines. Coast Eng 2014; 94(7): 10-22.
[31] Dodaro G, Tafarojnoruz A, Stefanucci F, et al. An experimental and numerical study on the spatial and temporal evolution of a scour hole downstream of a rigid bed. In: Proc. the International Conference on Fluvial Hydraulics, River Flow (2014); 3-5 September 2014; Lausanne, Switzerland. 2014; pp. 1415-22.
[32] Dodaro G, Tafarojnoruz A, Sciortino G, Adduce C, Calomino F, Gaudio R. Modified Einstein sediment transport method to simulate the local Scour evolution downstream of a rigid bed. J Hydraul Eng 2016; 142: 1-11.
[33] Ye Z, Wem D. Hydraulics. Beijing: China Communications Press 2001. (in Chinese)
[34] Ferziger J, Peric M. Computational method for fluid dynamics. 3rd ed. Berlin: Springer 2002.
[35] Melville W. Local Scour at Bridge Sites. Auckland: The University of Auckland 1975.
[36] Dey S, Barbhuiya A. Turbulent flow field in a scour hole at a semicircular abutment. Can J Civ Eng 2005; 32: 213-32.
[37] Calomino F, Tafarojnoruz A, De Marchis M, Gaudio R, Napoli E. Experimental and numerical study on the flow field and friction factor in a pressurized corrugated pipe. J Hydraul Eng 2015; 141(11)
[38] Gogus M. Effect of spur dike length on the horseshoe vortex system and the bed shear stress distribution. J Hydraul Res 2015; 53(2): 196-206.