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

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:

(1)

(2)

in which, x_i (i=1, 2, and 3) is the coordinate along axis in Cartesian coordinate system, u_i is the instantaneous velocity component in x_i direction; ρ is the density of clear water; t is time; u_jp is the velocity of the movement of the computational mesh; ν is the kinematic viscosity; τ_ij is the Reynolds stress defined as τ_ij = ν_t(∂u_i/∂x_j + ∂u_j/∂x_i)−(2/3)kδ_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 x_i 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:

(3)

(4)

where represents the generation of turbulence kinetic energy due to the mean velocity gradients; μ_t = ρC_μk² / ε; C_1ε, C_2ε, C_μ, σ_k, σ_ε are the empirical constants, with value of C_1ε = 1.44, C_2ε = 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:

(5)

In which subscript “p” stands for the center of control volume adjacent to wall; U_p 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:

(6)

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:

(7)

In which k_s⁺ = u*k_s / ν 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, k_s should be equal to the median grain size, d₅₀. At present, a value of k_s = d₅₀ 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:

(8)

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 U² / (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:

(9)

where γ is the sediments porosity; h is the elevation of the channel bed; and q_bi represents the bed load flux vector in the x_i 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 q_bi following the study by Qian and Wu[1]:

(10)

in which τ_i is the bed shear stress in x_i 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:

(11)

(12)

in which d_iis 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_* = d₅₀ [(s − 1)g / ν²]^1/3 is the dimensionless particle parameter; s = ρ_s/ ρ is the relative density; T is the non-dimensional excess bed shear stress; and T_m is the modified non-dimensional excess bed shear stress; T and T_m are expressed as:

(13)

(14)

where λ = e^{(0.45α + 0.2β)} is correction factor; α = 1 − σ_g; β = d₅₀ / d₉₀ − d₁₀ / d₅₀; σ_g = (d₈₄ / d₅₀ + d₅₀ / d₁₆)/ 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:

(15)

where φ is the angle of repose for sediments; C_L donates the lift coefficients for the particle; C_D is the drag coefficients for the particles; C_L/C_D=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:

(16)

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.

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

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.01712m³ 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 m³ 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.

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

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. (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 80^o 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. (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.

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.

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.

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, d_m, 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.

In addition, flows on two horizontal slices with z = -0.5d_m and z=0 (d_m: 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.5d_m 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.

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