## Abstract

Flow and scour around a vertical cylinder exposed to current are investigated by using a three-dimensional numerical model based on incompressible Reynolds-averaged Navier–Stokes equations. The model incorporates (i) *k*-*ω* turbulence closure, (ii) vortex-shedding processes, (iii) sediment transport (both bed and suspended load), as well as (iv) bed morphology. The influence of vortex shedding and suspended load on the scour are specifically investigated. For the selected geometry and flow conditions, it is found that the equilibrium scour depth is decreased by 50% when the suspended sediment transport is not accounted for. Alternatively, the effects of vortex shedding are found to be limited to the very early stage of the scour process. Flow features such as the horseshoe vortex, as well as lee-wake vortices, including their vertical frequency variation, are discussed. Large-scale counter-rotating streamwise phase-averaged vortices in the lee wake are likewise demonstrated via numerical flow visualization. These features are linked to scour around a vertical pile in a steady current.

## 1. Introduction

Flow and scour around a vertical cylinder have been investigated extensively over the past four decades, with specific relevance to the flow and scour around bridge piers in river engineering, and around pile foundations in marine and coastal engineering. While much has been written on this subject, comparatively few studies have been presented involving detailed three-dimensional numerical modelling of such processes. Important contributions to the numerical investigation of both flow and scour around vertical piles include [1–15].

The research presented in this paper is largely motivated by the previous work of Roulund *et al*. [3]. Their model was based on cohesionless sediment and simulated scour in the live-bed regime. However, unsteady flow structures such as lee-wake vortices (due to mesh resolution), as well as suspended sediment transport, were not incorporated in their simulations. Recently, Stahlmann [14] has investigated the effects of these processes on the scour case studied by Roulund *et al*. [3]. Their study further improved the results of scour depths achieved by Roulund *et al*. [3].

This study focuses on incorporating the effects of lee-wake vortices, as well as suspended sediment transport, within the three-dimensional numerical modelling of scour around a vertical pile. The fully coupled model provides a detailed overview of both phase-averaged (time-averaged over a certain number of periods of vortex shedding until no significant changes are observed in the phase-averaged quantities), as well as phase-resolved (instantaneous), flow structures in the lee wake of a vertical pile placed on a plane bed. Kirkil *et al*. [16] investigated coherent structures present in the flow field around a circular cylinder located in a scour hole with large-eddy simulation and laboratory-flume visualizations. Recently, Petersen [17] experimentally investigated the edge scour around scour protections of offshore wind turbine foundations, with special emphasis on the flow features in the lee-wake region. In their 2008 study, Kirkil *et al*. [16] addressed the presence of the phase-averaged streamwise vorticity fields in the near-wake region of the pile for a scoured bed. Petersen [17] emphasized the prominent effects of these vortices on the lee-wake scour around a scour protection of a vertical pile. Petersen showed that the edge scour at the downstream side of the scour protection is mainly caused by these vortices, two large-scale counter-rotating streamwise vortices in the lee wake, which have been visualized by sand. The effects of roughness (due to stone protection around the pile) on the formation of these vortices could not be interpreted, however, as a comparable study has yet to be undertaken with a rigid plane bed. Regarding the lee-wake vortices, Jacobsen *et al*. [15] have studied transient effects from vortex shedding on scour patterns. They have shown that the process of vortex shedding is affected by the sheared flow itself and the scour depth around the cylinder, and the shedding frequencies are found to vary along the height of the cylinder for various scour depths. In the present work, such a study will be directly carried out through three-dimensional numerical simulations of the flow around a cylinder placed on a rigid and plane smooth bed.

This paper is organized as follows. The hydrodynamic and turbulence models used are described in §2, and the sediment transport and morphological models are described in §3. Simulations regarding the flow features around a vertical circular pile placed on a rigid bed under steady flow conditions are described and discussed in §4. These will serve to validate the model as well as to help visualize phase-averaged flow features for the given flow conditions. Simulations of scour around a circular pile, to investigate the effect of unsteady flow structures and the suspended sediment transport, are subsequently presented and discussed in §5. Finally, conclusions are drawn in §6.

## 2. Hydrodynamic and turbulence models

### (a) Governing equations

In this section, a description of the computational model used throughout the present work is provided. The numerical model solves the incompressible Reynolds-averaged Navier–Stokes equations 2.1 where the mean strain-rate tensor is 2.2

These are combined with the local continuity equation
2.3
Here, *u*_{i} are the mean (phase-resolved) velocities, *x*_{i} are the Cartesian coordinates, *t* is time, *p* is the pressure, *ν* is the fluid kinematic viscosity, *ρ* is the fluid density and *τ*_{ij} is the Reynolds stress tensor, which accounts for additional (normal and shear) stresses due to momentum transfer from turbulent fluctuations. Throughout the present work, the Reynolds stress tensor will be defined according to the constitutive relation
2.4
where *δ*_{ij} is the Kronecker delta, *ν*_{T} is the eddy viscosity,
2.5
is the turbulent kinetic energy density, and the overbar denotes time averaging. To achieve closure, the two-equation *k*-*ω* turbulence model of [18,19] is adopted. In this model, the eddy viscosity is defined by
2.6
which incorporates a stress-limiting feature, with . This model additionally uses transport equations for the turbulent kinetic energy density *k*
2.7
as well as for the specific dissipation rate *ω*
where
2.8
and is the Heaviside step function, taking a value of zero when the argument is negative and a value of unity otherwise. The standard model closure coefficients are used: , , , , and , *β*=*β*_{0}*f*_{β} and *β*_{0}=0.0708, with
2.9

### (b) Boundary conditions

The hydrodynamic model described above is subject to the following boundary conditions. At friction wall boundaries, a no-slip condition is imposed whereby velocities are set to zero. Alternatively, Neumann conditions are applied to the three components of the velocity and scalar hydrodynamic quantities at the symmetry boundaries, which is in line with the approach in [3]. The ‘lid’ at the top boundary is not rigid and finite velocities are allowed perpendicular to the upper boundary. The application of a ‘lid’ at the top boundary rather than a free surface facility is a commonly used and valid approach for the flow conditions without large free surface disturbances as in the case of wave modelling cases. Regarding the effect of the ‘lid’, this issue has been discussed in [3] in conjunction with their numerical simulations of flow and scour around a pile. With a free surface, and for sufficiently large Froude numbers, *Fr*=*U*/(*gh*)^{1/2}, in which *h* is the height of the computational domain, *U* is the mean (undisturbed depth-averaged) flow velocity and *g* is the gravitational acceleration, the free surface will, in real life, exhibit some variation in the vicinity of the cylinder; there will be a run-up in front of the cylinder and a depression around the side edge and at the back of the cylinder. Roulund *et al*. [3], pp. 373–378 pointed out that this variation ceases to exist for Froude numbers smaller than *O*(0.2). Now, in the present case, the Froude number, *Fr*=*O*(0.4), is not radically different from *O*(0.2), and therefore the flow will be fairly well represented in the lid simulations although some variations may be expected such as the absence of a flow in the radial direction, caused, otherwise, by the head difference between the surface elevation in front and at the side edge of the cylinder. Regarding the side boundaries, for all the quantities except the velocity, *U*, zero-gradient boundary conditions are applied, and for the velocity the slip type boundary condition, in which the normal component of the velocity is taken to be zero, is applied.

To drive the flow conditions within the model, prior to the actual simulation, for the same calculation domain without the cylinder, the flow is first driven by a horizontal pressure gradient specified as a constant based on the desired undisturbed friction velocity, *U*_{f}=[−(*h*/*ρ*)×(∂*p*/∂*x*)]^{0.5}. As the steady-state flow condition is reached, the *U*, *k* and *ω* profiles are taken to be used as the inlet boundary conditions in the actual simulation for the calculation domain with the cylinder. At the inlet, for the remaining quantities, Neumann boundary conditions are applied. At the opposite right-hand boundary (outlet), zero-gradient boundary conditions are imposed for all quantities except the pressure term, *p*, which has been taken as zero at the outlet boundary.

At the bottom as well as the cylinder, a generalized wall function approach developed by Fuhrman *et al*. [20] is used. Under this approach, the friction velocity *U*_{f} is determined from the tangential velocity at the nearest cell centre based on the profile of Cebeci & Chang [21],
2.10
2.11
2.12
Cebeci & Chang [21] generalized the van Driest [22] profile to incorporate potential roughness effects. Here, *y*_{c}=*Δy*/2 is the normal distance from the wall to the cell centre, expressed in terms of wall coordinates as , *Δy* is the near-wall cell thickness, *κ*=0.4 is the von Karman constant and *k*_{s} is Nikuradse's equivalent sand grain roughness. At the bed *k*_{s}=2.5*d* is used, whereas at the cylinder *k*_{s} is kept sufficiently small (*k*_{s}=1×10^{−5} m) such that it is effectively modelled as hydraulically smooth. Here, *d* is the grain size. The friction velocity is then used within the generalized wall functions for *k* and *ω* in the cells nearest to the wall, according to
2.13
and
2.14
The first arguments in these functions ensure that these variables follow their proper scaling *k*∼*y*^{2} and *ω*∼1/*y*^{2} for near-wall cells within the viscous sub-layer (e.g. [18]). The values and are used, which ensure a continuous transition to the (fully turbulent) second arguments at *Δy*^{+}=*δ*^{+}, where *δ*^{+}=11.626 is taken as the viscous sub-layer thickness (in dimensionless wall coordinates).

### (c) Computational mesh

The calculation domain is discretized into finite volumes of quadrilateral blocks in varying shapes and dimensions. Figure 1 shows the boundaries of the domain on an example of computational mesh used for scour calculations. During the live-bed calculations, the mesh is continuously updated to adjust the changes of bed topography. For further details, see §3*b*.

The calculation domain has the following dimensions: length, *l*=20*D*, width, *w*=15*D*, and height, *h*=2*D* (or 3*D*), in which *D* is the pile diameter. The total number of cells is *O*(10^{5}). Table 1 summarizes the mesh characteristics for different domain heights and calculation types. The pile is located at the centre of the domain (*x*=0 and *y*=0). As an indication of computational expense, a fully coupled hydrodynamic and morphological calculation lasting 1 min of physical time for the three-dimensional scour around a vertical pile problem requires approximately 10 days of CPU time, when simulated in parallel on eight modern processors. To give an indication of the CPU times, the authors would like to note that, for the same computational mesh and using the same computational power, only the hydrodynamic calculations for 1 min of physical time last approximately 1.5 days, and 3 days if the suspended sediment calculations are also included. In all the calculations, the time step is kept as variable to ensure that certain Courant numbers selected separately for the numerical stability in each of hydrodynamic, suspended sediment and morphological calculations are not exceeded. Throughout the study, it has been seen that the morphological calculations have governed the speed of the simulations. The authors also would like to note that considerable effort has been put into optimizing the computational mesh for the best grid and time convergence. For the scour simulations, finer computational meshes, as used in rigid bed simulation with *h*=3*D*, could not be used as CPU times increased significantly due to numerical stability considerations.

In this study, basically, two kinds of calculations are carried out: the rigid bed calculations to investigate the flow features only; and scour calculations. In rigid bed calculations, the bottom of the calculation domain is kept rigid by turning off the morphological model, whereas in scour calculations the bed is continuously updated. The rigid bed simulations are carried out for two different domain heights (*h*=2*D* and 3*D*), whereas in the scour simulations only the calculation mesh with *h*=2*D* is used.

## 3. Sediment transport and morphological models

This section describes the sediment transport and morphological models used herein. As the implementation of these models, including a detailed account of numerical aspects, has been described in detail by Jacobsen [23], as well as in the recent publication of Jacobsen *et al*. [24], only essential details are provided in what follows.

### (a) Sediment transport model

The model makes use of separate bed and suspended sediment load descriptions. The rate of bed load sediment transport *q*_{B} is calculated based on the method described in detail by Engelund & Fredsoe [25], as well as Roulund *et al.* [3], who generalized the well-known transport formulation of Engelund & Fredsoe [26] to account for three-dimensional effects as well as bed-slope modifications.

The suspended sediment model is, alternatively, based on a turbulent-diffusion equation (e.g. [27], p. 238) of the form
3.1
where *c* is the suspended sediment concentration, *w*_{s} is the settling velocity, and *β*_{s}=1 is used throughout, i.e. the sediment diffusivity is taken as equal to the eddy viscosity. The settling velocity is solved for a given grain diameter *d* using the drag-coefficient approach described in, for example, [27]. Equation (3.1) is solved on a sub-set of the main computational mesh, i.e. where near-bed cells below a so-called reference level *b* are removed, as described in detail by Jacobsen *et al*. [24], fig. 2.

In all boundaries except the bottom boundary, a zero-flux condition for *c* is used. At the bottom boundary so-called reference concentration boundary conditions are imposed. More specifically, the method of Engelund & Fredsoe [26] is used, with the concentration at the reference level set to
3.2
where *c*_{0}=0.6, and the linear concentration *λ*_{b} is
3.3
where
3.4
is the probability of moving grains, and
3.5
is the Shields parameter, in which *τ*_{0}=*ρU*_{f}^{2} is the wall shear stress, *s*=*ρ*_{s}/*ρ* is the specific gravity of the sediment grains, and *g* is the gravitational acceleration. Throughout the present work, the coefficient of dynamic friction is set to *μ*_{d}=0.51. The critical Shields parameter *θ*_{c} is computed from a base value *θ*_{c0}=0.05, accounting for bed-slope effects as in [3]. Following [20], throughout the present work, the reference level is taken as *b*=*α*_{1}*d*=3.5*d*.

To prevent un-physical ‘overloading’ conditions (i.e. where reference *c*_{b} is forced to be smaller than the concentration immediately above) from occurring in the model, the solution suggested by Justesen *et al*. [28] is finally invoked. That is, if the concentration close to the bed exceeds the reference concentration *c*_{b}(*θ*), the value in practice is taken from the cell adjacent to the boundary.

### (b) Morphological model

The morphology of the bed elevation *h*_{b} is based on the sediment continuity (Exner) equation
3.6
where *n*=0.4 is the bed porosity, and D and E are the deposition and erosion stemming from the suspended sediment model, respectively. Further details of the evaluation of deposition and erosion terms in the above equation are given in [24], including corrections for the non-sloping beds. It is stressed that the simulated bed morphology within the present work is continuous, always being based on the instantaneous sediment transport fields, i.e. the model does not make use of morphological rates averaged over a certain duration or any other time scale. Accordingly, morphological and hydrodynamic times are equivalent. The temporal integration of (3.6) is herein based on the explicit Euler method. For more specific details regarding the numerical evaluation of the three terms within the right-hand side of (3.6), the interested reader is referred to [24]. It is finally noted that, to prevent excessive erosion induced by the imposed uniform flow at the boundary, the sea bed is fixed at the left (inlet) boundary and relaxed to full morphology over a distance spanning a few pile diameters.

Experience has shown that, if left unchecked, the morphological model can lead to local bed slopes in excess of the angle of repose. To combat such un-physical steepening, the physically based sand slide model described in detail in [3] is implemented. In the present work, this sand slide model is activated at positions where the local bed angle exceeds the angle of repose *ϕ*_{s}=32° and is de-activated once the local bed angle is reduced to 30.0°. Additionally, some local filtering of the bed was necessary in the scour runs for stability reasons. For this purpose, a filtering algorithm based on the least-squares method is applied to the morphological rates in the near vicinity of the pile to smooth out small-sized (i.e. high wavenumber) ripples that can occur inside the scour hole adjacent to the pile, where the numerical stability of the mesh motion computations becomes critical due to increased rates of change in bed elevations. The strength of filtering inside and around the scour hole is kept variable with the use of a hyperbolic function based on the radial distance from the origin of the domain, so that the far-field bed ripples are resolved but the small-sized ripples inside the scour hole are filtered.

The equations constituting the fully coupled model outlined above are solved numerically using the open-source CFD toolbox OpenFOAM, v. 1.6-ext, making use of a finite volume spatial discretization with a collocated variable arrangement, in conjunction with a standard PIMPLE algorithm. Again, for further details see [24]. For the simulation in which the vortex shedding is turned off, more diffusive schemes are used for the divergence terms in (2.1).

## 4. Rigid bed simulations

In this section, the flow features around a vertical cylindrical pile mounted on a horizontal plane soil bed and subjected to a steady current are investigated through numerical simulations for the two meshes listed in table 1. For the simulations in this section, the bed is kept rigid by turning off the morphological model. The test conditions of the simulations (for both rigid bed and scour simulations) with respect to domain height are given in table 2. As shown previously in table 1, the rigid bed simulations are carried out for two different domain heights (*h*=2*D* and 3*D*), whereas in the scour simulations only the computational mesh with *h*=2*D* is used.

### (a) Horseshoe vortex

The main flow features around a vertical cylindrical pile exposed to a steady current on a plane bed are the horseshoe vortex formed in front of the pile, the lee-side vortices (usually in the form of vortex shedding), contracted streamlines around the sides of the pile and the down-flow in front of the pile due to deceleration of the flow. The horseshoe vortex in front of the pile is an end result of the adverse pressure gradient on the bed upstream of the pile, which is actually induced by the presence of the pile in a steady current. This adverse pressure gradient enforces a three-dimensional boundary layer separation just in front of the pile and results in a vortical structure surrounding the upstream perimeter of the pile and fading away downstream. The size of the horseshoe vortex in front of a circular pile with a smooth bed depends mainly on [29] the ratio of the bed-boundary-layer thickness to the pile diameter, *δ*/*D*; the pile Reynolds number, *Re*_{D}=*UD*/*ν* (or alternatively the bed-boundary-layer-thickness Reynolds number, *Re*_{δ}=*Uδ*/*ν*), where *U* is the free-stream velocity; and the ratio of bed roughness to pile diameter, *k*_{s}/*D*, in the case of rough beds.

Figure 2 compares the amplification of the bed shear stress (*α*_{τ}) along the *x*-axis calculated from the model (for *h*=2*D*) with that measured in the smooth bed experiments by Roulund *et al*. [3] and Hjorth [30] (see fig. 6.18.b on p. 133 of [30] and the Test-1 data in fig. 16 on p. 376 of [3]). Here, the amplification of the bed shear stress is defined as , in which *τ*_{0} is the bed shear stress and is the undisturbed bed shear stress. The negative bed shear stress corresponds to the location of the horseshoe vortex in front of the pile. Figure 2 shows that the numerical and experimental results agree well outside the horseshoe vortex for *x*/*D*<−1.2. However, for −1.2<*x*/*D*<−0.5, the bed shear stress is apparently under-predicted by the numerical model with respect to the data given by Hjorth [30], the maximum difference between the model and the experiment being around 30%. The latter is consistent with the findings given in [3]. No clear explanation has been given in [3] or found for the observed discrepancy between the numerical model and the experiments.

Figure 3 compares the amplification of the bed shear stress defined as the ratio of magnitude of the bed shear stress vector to the value of the undisturbed bed shear stress () obtained from the present model (for *h*=2*D*) and that from [30] smooth-bed experiments (see fig. 6.18.b on p. 133 of [30]). The present calculations and the measurements are seen to be in good agreement in terms of the magnitude of the bed shear stresses around the cylinder. However, the radial location of the maximum bed shear stress is found to be rather different from the measurements, which might be a possible explanation for the difference between the model and Hjorth's results [30] in figure 2.

### (b) Lee-wake flow

The lee-wake flow around the cylinder depends mainly on the pile Reynolds number, *Re*_{D}. For *Re*_{D}>40 (e.g. [31]), the flow in the lee wake becomes unsteady and the vortex-shedding regime begins to occur with vortices shed alternately at either side of the pile at a certain frequency. In this study, *Re*_{D} is equal to 1.7×10^{4}, for which the wake flow is given as completely turbulent by Sumer & Fredsoe [31]. These vortices are caused mainly by the rotation in the boundary layer over the surface of the pile. The shear layers emanating from the side edges of the pile roll up to form these vortices in the lee wake of the pile. The emanating position of these shear layers on the cylinder, the so-called position of boundary layer separation, changes over the height of the cylinder when the velocity profile is not uniform over the depth of the domain.

To investigate the flow features in the lee wake, first, the position of the boundary layer separation over the height of the cylinder is drawn for the rigid bed numerical model simulations in which the bed roughness is in transitional range (close to hydraulically smooth conditions) and the cylinder roughness is taken as hydraulically smooth (table 2). Figure 4 shows the side view of the position of boundary layer separation around the cylinder for both computational meshes (*h*=2*D* and 3*D*), where the velocity on the cylinder surface becomes zero or negative with respect to the streamwise direction. It can be seen in figure 4 that the boundary layer separation close to the bottom starts earlier as the bottom velocities are smaller in magnitude than the upper section of the domain, and the bed roughness is almost smooth and does not result in a delay in boundary layer separation close to the bottom as in the case of flow around a vertical hydraulically smooth cylinder placed on a rough bottom. As in the latter case, the rough bottom increases the exchange of momentum between the outer flow region and the wall, and, thus, in the streamwise momentum in the boundary layer resulting in a delay in the separation [32]. As the velocities close to the bottom are smaller in magnitude, the magnitude of the vorticity supply in the near lee-wake region becomes less than the upper parts of the domain. As an important implication of this, the vortex-shedding mechanism slows down close to the bottom, causing a delay in the shedding regime at the bottom, and, thus, the shed vortical structures along the height of the domain are broken apart into two cells in the vertical direction.

Figures 5 and 6 give sequences of power spectra of the lift force applied on the cylinder in the *y*-direction obtained for different values of *z*/*h*, for the computational mesh with *h*=2*D* and 3*D*, respectively. Here, *ϕ*_{L} and *σ*_{L}^{2} are the power spectrum and the variance of the lift force, respectively, and *S*_{t} is the Strouhal number (*S*_{t}=*fU*/*D*), where *f* is the vortex-shedding frequency, *U* is the mean flow velocity and *D* is the pile diameter. Both simulations show a clear shift from lower to higher *S*_{t} as *z*/*h* increases, indicating more frequent shedding of vortices further from the bed. The shift in Strouhal number clearly shows that the vortex shedding occurs in two separate cells, rather than one cell of domain height. This shift is observed to be around *z*/*h*=*O*(0.2) for both domains (*h*=2*D* and 3*D*). Regarding the effect of domain height on the scour, Sumer & Fredsoe [33], fig. 3.26 show that the boundary layer thickness has a clear effect on the scour depth up to a boundary layer thickness to pile diameter ratio equal to 4. Also, it can be seen from figure 5 that the spectral peaks for *h*=2*D* are rather irregular (composed of multiple peaks rather than just one peak at one frequency) compared with the ones given for *h*=3*D*. However, the highlight of this section might be given as the behaviour of vortex-shedding changes close to the bottom due to the bottom shear effect. One final note with regard to the vortex-shedding frequency is that the values of the Strouhal number are apparently different from the familiar value 0.2. This is linked with the fact that the mean flow velocity *U* is used when normalizing the frequency, rather than the ‘local’ approach velocity.

The two-cell structures of the shed vortices are much more visible in figure 7, in which the magnitude of the vorticity in the vertical direction (*ω*_{z}=∂*u*_{y}/∂*x*−∂*u*_{x}/∂*y*, in which *ω*_{z} is the vorticity in the vertical direction, *u*_{x} is the velocity component in the *x*-direction and *u*_{y} is the velocity component in the *y*-direction) around the cylinder at a selected instant in time is given (for *h*=3*D*).

The upper and lower bounds of the colour scale of figure 7 are limited to small values of vorticity to make the opposite fields easily distinguishable. As seen from figure 7, the vertical structure of the vorticity in the *z*-direction is broken simply into two parts close to the bottom, and the shed vortical structure looks like either a tornado with a tail that follows the main body from behind or a broken tornado where the top and bottom parts rotate in opposite directions with respect to each other.

The flow around the cylinder over a rigid bed can be further visualized in an accompanying video (see electronic supplementary material). The video is composed of four views, each presenting the motion of a different quantity to visualize the coherent structures in the flow. In the upper left and right quadrants of the video, the Q-criterion [34] is used to educe the main coherent structures such as the horseshoe vortex and the lee-wake vortices. Here, an overall look at the vortical structures around the cylinder is given. One can see the horseshoe vortex formed in front of the pile and the lee-wake vortices that are dispersed in the flow as they progress. In the upper right view, the threshold for *Q* is set to higher values so that a detailed picture of lee-wake vortices just behind the cylinder is given. Here, it is easily seen that the vortex shedding occurs in two cells, where the bottom cell extends approximately up to *z*/*D*=1.0 from the bottom. In the lower left view of the video, iso-surfaces of vorticity in the *z*-direction are visualized. In this view, it can be seen that the bottom part of the vortical structure moves more slowly than the upper part, like the swing of a skirt, and that the shedding occurs more slowly close to the bottom. The lower right view shows iso-surfaces of vorticity in the streamwise direction. Here, again, the two-cell structure of vortex shedding is prominent and it gives hints about the formation of large-scale counter-rotating streamwise phase-averaged vortices (LSCSVs), which are discussed further in detail below.

Figure 8 depicts a sequence of bed shear stress (*τ*_{0}) vectors and contour lines of magnitude of bed shear stress amplification () at the bottom (for *h*=2*D*), obtained from the unsteady solution of the present model, over one-half period of vortex shedding. Here, the unsteady behaviour of the wake is evident. Incidentally, the figure further indicates the presence of the horseshoe vortex.

The present model compares well with the pressure measurements of Dargahi [35] as given in [3]; these particular results are not repeated here for brevity. Figure 9 shows the phase-averaged pressure field around the cylinder along the symmetry plane (for *h*=3*D*). The quantity depicted is the ratio of the pressure coefficient, , at any point in the domain, to the pressure coefficient computed at the stagnation point of the top boundary, (at *x*=−0.5*D*, *y*=0 and *z*=*h*), in which *p* is the pressure (, where is the time-averaged pressure and *p*_{0} is the hydrostatic pressure) and *U* is the mean flow velocity. The colour scale in the figure is selected to make the pressure gradients clearly visible. One can see from this figure that there exists a pressure gradient in the vertical direction over the symmetry plane. This pressure gradient drives a mean flow in the upward direction. The same would hold true if the phase-resolved pressure values are used. This is simply due to the far-field velocity profile which creates a higher pressure field close to the bottom. Similarly, there exists a pressure gradient at the top boundary driving again a mean flow towards the side boundaries at the relatively far field with respect to the cylinder.

Figure 10 depicts the pressure field around the cylinder (for *h*=3*D*) at the horizontal cross section corresponding to *z*=1*D*. The quantity depicted is the same as that in figure 9. One can again see a pressure gradient which drives a mean flow towards the centre of the domain in the transverse direction. The mean pressure gradients in the relatively far-field of the cylinder presented in figures 9 and 10 creates LSCSVs downstream of the cylinder.

Figure 11 shows phase-averaged velocity streamlines over transverse planes at *x*=0.75*D*, *x*=2.5*D*, *x*=5.5*D* and *x*=8.5*D* (for *h*=3*D*). As seen, LSCSVs become much more evident as our focus moves further away from the cylinder. At *x*=0.75*D*, one can see a region in the middle with upward directed velocity streamlines. This is the region with a strong vertical pressure gradient just behind the cylinder (figure 9). As the shed vortices disperse and the velocity field downstream becomes more uniform, secondary circulation cells disappear, while the two large circulation cells remain.

Figure 12 depicts the LSCSVs, visualizing the phase-averaged velocity streamlines in figure 12*a*, and the iso-surfaces of the phase-averaged vorticity in the streamwise direction (*x*) in figure 12*b* (for *h*=3*D*). To visualize the LSCSVs in the form of streamlines in figure 12*a*, the phase-averaged velocities are used. As the flow is mainly unidirectional and the velocity component in the streamwise (*x*) direction is much stronger than the other two components (in the *y* and *z* directions) in most parts of the domain, here, the phase-averaged velocity component in the *x*-direction only is reduced by multiplying by a factor of 0.02, which is found by trial and error, and the other two components are kept as they are. By doing so, the phase-averaged flow is slowed down in the streamwise direction only, so that the LSCSVs can be presented in the form of spirals in the three-dimensional domain. In addition to the LSCSVs shown in figure 12, another pair of vortices close to the bottom and adjacent to the cylinder rotating in the opposite direction merits attention. These vortices are partly due to the horseshoe vortex and the change in the vortex-shedding frequency at the bottom, resulting in a delay of almost one shedding period of the upper cell.

It should be noted that, as mentioned previously, recently Petersen [17] studied edge scour at the scour protection around the foundation of a wind turbine by physical model experiments. He found that the edge scour at the downstream side of the scour protection is mainly caused by the LSCSVs in the lee wake of scour protection and the pile system. These vortices were also clearly visualized by his video taken downstream where the LSCSVs were visualized by the sand. The bed morphology obtained after the experiments in [17] consisted of two elongated scour holes with a longitudinal ridge between them, and the bed form (associated with the edge scour) appears to be the ‘signature’ of the LSCSVs.

## 5. Simulation of the scour process

Scour around a vertical pile structure has been previously investigated in detail in [3], where steady-state flow calculations around a vertical pile over an initially plane bed were incorporated with a morphology model based on bed load sediment transport and sand slide processes only. In this study, the additional effects of unsteady flow structures and suspended sediment transport are investigated for the test conditions (for *h*=2*D*) indicated in table 2. For this purpose, three different scenarios with the same flow conditions and the same computational mesh are studied. In scenario 1, both the unsteady flow features (vortex shedding) and the suspended sediment transport are turned on. In scenario 2, only the vortex shedding is turned on and the suspended sediment transport is turned off. In scenario 3, both the vortex shedding and suspended sediment transport are turned off. The three scenarios are summarized in table 3. Note again that vortex shedding can be conveniently and effectively switched off by using more diffusive schemes (e.g. Gauss upwind) for the divergence terms in (2.1). The computational mesh used in the scour simulations is with the domain height, *h*=2*D* (table 1). The flow conditions in the simulations are given in table 2 for *h*=2*D*.

The undisturbed Shields parameter in the test conditions is given as *θ*=0.13, which is larger than the critical value for the initiation of motion, *θ*_{c}=0.05, implying that the scour is in the live-bed regime. The calculated maximum Shields parameter around the cylinder is around *θ*≈1.3; therefore, suspended sediment transport can be expected to be active in scour around the cylinder. In the simulations, bed ripples are resolved. No extra treatments are applied to trigger the formation of the bed ripples, i.e. they arise naturally within the computational domain as the ripples form as a result of the physical instability of the sediment bed when exposed to a steady boundary layer flow provided that the grain Reynolds number, d*U*_{f}/*ν*, is smaller than *O*(20) [36]. The sediment transport description in the present model is based on uniform-size sand, i.e. any effects of sediment gradation are not accounted for.

Figure 13 shows a sequence of pictures illustrating the time evolution of the scour hole obtained in scenario 1. In the figure, one can see (i) the semi-circular shape (in plan view) of the upstream part of the scour hole with a slope equal to the angle of repose, (ii) the formation of a ‘bar’ downstream of the pile (the deposited sand), and its downstream migration with the contribution of LSCSVs, (iii) the formation of a gentler slope of the downstream side of the scour hole, (iv) the formation and migration of bed ripples at the downstream side of the scour hole, and, finally, (v) the bed material eroded and piled up along the edges of the scour hole downstream by means of the counter-rotating streamwise phase-averaged vortices close to the bottom shown in figure 12*a*.

The equilibrium scour depth obtained in scenario 1 is *S*/*D*=0.91, where *S* is the scour depth in front of the pile. This value agrees well with the existing data, for the present water-depth-to-pile-diameter ratio *h*/*D*=2.0 (see fig. 3.26, p. 180, in [33]) in the simulation and the ratio between the approach velocity and the critical velocity for initiation of sediment motion is *U*/*U*_{cr}≈1.6 (see fig. 3.25, p. 179, in [33]).

Figure 14 shows the time development of the scour depth at the upstream side (figure 14*a*) and at the downstream side (figure 14*b*) of the pile for the three scenarios. Figure 15, on the other hand, shows the scour-hole profiles along the upstream–downstream symmetry line after 2 min of simulation time, again for all the three scenarios. In figure 14*a*, one can see that the upstream equilibrium scour depth decreases by almost 50% when the suspended sediment transport is not incorporated in the computations, regardless of vortex shedding being enabled or disabled. In figure 14*b*, the equilibrium scour depth ratio at the downstream side of the pile reaches to *S*/*D*=0.52 for scenario 1. For scenarios 2 and 3, however, the scour depth ratio reaches to *S*/*D*=−0.23 and *S*/*D*=−0.16, respectively, the negative values indicating accretion rather than scour. This is because the sediment could not be picked up and carried further away by the flow, but instead piled up adjacent to the pile downstream. The initiation of scour at the downstream side of the pile for scenario 2 within the first 10 s might be attributed to the larger bed shear stresses under the shed vortices.

In figure 15, the slope of the scour hole at the upstream side of the pile is equal to the internal angle of friction of the sediment grains, and the downstream slope is relatively less steep, which agrees well with the observations [33]. The formation of the ‘bar’ accumulated adjacent to the pile in scenarios 2 and 3 has been explained previously. The ‘bar’ formed at the downstream side in scenario 1, which is further away than scenarios 2 and 3, is in agreement with the observations of scour around a vertical pile in steady current.

The present simulation results for scenario 3 differ from the results in [3], where the simulation conditions resemble scenario 3. Two major differences may be noted, namely (i) the equilibrium scour is, in the present simulation, reached much faster than in [3] and (ii) the equilibrium scour depth is a factor of 2 smaller than in [3] although the Shields parameter and the boundary-layer-thickness-to-diameter ratio of the present simulation are much the same as in [3]. No clear explanation has been found for this discrepancy. Nevertheless, this is not a problem as the main objective of the present scour simulations is to investigate the influence of the suspended load as well as the vortex shedding on the scour process, and the present simulations serve the purpose.

At this juncture, it is also interesting to elaborate on the influence of the suspended load on scour. Laboratory observations in the previous research indicated that scour around a pile in steady currents is not influenced very significantly by the presence of suspended load (e.g. fig. 3.25 in [33]). The latter source shows that the scour depth increases by a factor of 1.3 (for a nearly uniform size sediment) when the Shields parameter is increased to a value corresponding to the suspension regime including the sheet flow transport. By contrast, the present simulations indicate that the increase in the scour depth is a factor of 2 when the suspended load is switched on. This aspect of the problem appears to be worth exploring in future research; see also the discussion in [37].

## 6. Conclusion

A three-dimensional hydrodynamic model, incorporated with a *k*-*ω* turbulence model, sediment transport (bed and suspended load) description and a morphological model is used to investigate the flow and scour around a vertical circular pile exposed to a steady current. The model, tested and validated against the experimental data in the literature, has been used to study flow processes to visualize the phase-averaged flow structures and to study the influence of the unsteady flow features and incorporation of the suspended sediment transport on the scour around a vertical pile.

The numerical model results show that the vortex shedding around a vertical pile placed on a plane bed occurs namely in a two-cell fashion, one cell very close to the bottom *z*<1.0*D* and the other one above the first cell, of which the shedding frequency is almost doubled with respect to the bottom. The variation of the shedding frequency over the height of the pile is found to be due to the early boundary layer separation on the surface of the pile close to the bottom because of bottom roughness, a smooth wall and the decrease in the magnitude of the vorticity feeding the lee wake as a result of smaller velocities at the bottom.

The phase-averaged flow features were first studied by visualizing the phase-averaged pressure fields in the lee-side of the pile. Here, a vertical pressure gradient continuous along the symmetry plane normal to the transverse direction at *y*=0*D*, strongest at −0.5<*x*/*D*<−1.5 and driving an upward-directed mean flow, is observed. Additionally, there exist horizontal pressure gradients close to the bottom and top boundaries that drive horizontal mean flows towards the centre close to the bottom and towards the sides of the domain close to the top boundary. The formation of the LSCSVs in the far-wake region is attributed to these mean pressure gradients. Rather small-scale streamwise phase-averaged vortices adjacent to the pile close to the bottom at the near-wake region could be attributed to the inherited vorticities from the horseshoe vortex in front of the pile and the increase in the vortex-shedding period (almost doubling the period in the upper cell) close to the bottom. The presence of the phase-averaged streamwise vortices is visualized through streamlines and the iso-surface of the vorticity magnitude in the streamwise direction.

The influence of unsteady flow structures (vortex shedding) and suspended sediment transport on the scour around a vertical pile has been investigated by considering three scenarios, in which these individual features are systematically switched on and off. When using the full model, i.e. with both vortex shedding and suspended load enabled, an equilibrium scour depth upstream of the pile *S*/*D*=0.91 has been obtained, which agrees well with existing data [33]. Alternatively, the upstream equilibrium scour depth is decreased by almost 50% when the suspended sediment transport is not incorporated in the simulations, regardless of whether vortex shedding is enabled or disabled.

With the full model, it can be seen that the bed material is eroded and piled up along the edges of the scour hole downstream by means of the small-scale counter-rotating streamwise phase-averaged vortices close to the bed, while the LSCSVs contribute to the formation of a ‘bar’ downstream of the pile (the deposited sand) and its downstream migration. The slopes of the scour hole at both the upstream and downstream sides of the pile in the simulations are found to be in agreement with observations in [33].

Alternatively, turning off the suspended sediment transport decreased the overall transport and resulted in accumulation at the downstream perimeter of the pile for the first 120 s of simulation. The initiation of scour downstream of the pile within the first 10 s, where the vortex shedding is active but the suspended transport is not, is attributed to the larger bed shear stresses under the shed vortices. In the cases considered, including or excluding vortex shedding has had only a minimal impact within the duration of the simulations.

## Funding statement

C.B., B.M.S. and D.R.F. acknowledge support from FP7-ENV-2013.604-3 European Union project ASTARTE (Assessment, STrategy And Risk Reduction for Tsunamis in Europe), grant no. 603839. C.B. also acknowledges the support of a post-doctoral grant from The Scientific and Technical Research Council of Turkey (TUBITAK, grant no. 2219). C.B. and B.M.S. additionally acknowledge Innovative Multi-purpose Offshore Platforms: Planning, Design and Operation (MERMAID), 2012–2016, grant agreement no. 288710 of the European Commission 7th Framework Programme for Research.

## Footnotes

One contribution of 12 to a Theme Issue ‘Advances in fluid mechanics for offshore engineering: a modelling perspective’.

- © 2014 The Author(s) Published by the Royal Society. All rights reserved.