Three-Dimensional Rans Modeling of Flow Around Circular Piers using Nested Grids

Abstract The paper introduces a three-dimensional numerical model that solves the Reynolds Averaged Navier-Stokes (RANS) equations on a curvilinear grid system using a novel nested grid approach. The main benefit of the model is the possibility to model locally complex hydraulic features in large rivers like the flow field at hydraulic structures. The entire study domain in such a case can be discretized with a coarser resolution, whereas a much finer resolution can be applied to a defined zone of the obstructions, where a detailed description of the flow field is needed. The model is tested on a laboratory experiment carried out at the Georgia Institute of Technology, where the flow field around a single and two double circular cylinders in a flatbed flume was studied. Simulated flow velocity, turbulent kinetic energy and bed shear stress distributions are in good agreement with measurements. However, deviations downstream of the piers indicate the limitation of the steady state description of the flow in the unstable wake zone. Nevertheless, the nested grid approach presented herein is a promising step towards the modeling of the local scouring phenomenon due to the relatively low computational demand.


INTRODUCTION
Local scouring around hydraulic structures can lead to their failure, hence, knowledge of the mechanism of scouring and the prediction of the characteristic scour geometry are of great importance during the design of such structures. Investigation of the flow pattern around bridge piers, abutments, spur dikes etc., is, therefore, a topic of research in the field of river engineering. The flow around a pier is considerably complex, three-dimensional and turbulent. Characteristic components of the flow are (Graf and Istiarto, 2002): strong downflow at the face of the pier, the formation of the horseshoe vortex around the pier, which is caused by the influence of adverse pressure gradients on the approaching boundary layer, and upwelling vortex shedding downstream of the pier. The reproduction of this complex flow field with computational modeling is a challenging task for researchers and practitioners. Moreover, numerical modeling of the flow structure and its interaction with the erodible river bed requires high computational effort. Hence, until recently only physical model studies were carried out. Some of the most relevant laboratory experiments on pier scour were reported by Shen et al. (1969), Melville and Raudkivi (1977), Jain (1981), Raudkivi and Ettema (1983), Dargahi (1990), Dey (1997), Melville and Coleman (2000) and Oliveto and Hager (2002). Physical model investigations mostly result in formulas for the maximum scour depth and/or the temporal behavior of scouring. However, generalized formulas are based on idealized flow conditions, like steady flow in a rectangular channel with homogeneous bed material. Therefore, those formulas are not suitable in all river engineering cases with a more complex geometry. Moreover, the formulas often result in different solutions for the same flow situation. For unique problems the use of small-scale physical models can be an alternative. These models are time consuming and costly, though. Due to increasing computer performance, CFD models have become suitable tools to carry out flow analysis for different river engineering problems. However, their applicability has to be proven using laboratory or field measurements.
A key factor for successful morphodynamical simulations is the accurate estimation of the flow field and bed shear stress distribution. In this study a numerical model that solves the Reynolds Averaged Navier-Stokes (RANS) equations is used to study the flow and turbulence field around bridge piers. A novel nested grid algorithm is implemented in the RANS model to calculate the hydrodynamics in a cost-effective manner. Model validation is carried out using experimental data. Three flume experiments that were conducted at the Georgia Institute of Technology are chosen. The experiments are flow around a single pier (case A) flow around two piers that were arranged in different patterns (case B and C). Velocity measurements collected through an Acoustic Doppler Velocimeter (ADV) are available for comparison. For cases B and C the two piers were arranged 4 pier diameters apart (measured center to center) at angles of 30 and 45 degrees to investigate how their alignment affects the flow (Fig. 1).

Fig.1
Layout of studied pier arrangements.

EARLIER WORK
This paper deals with flow and turbulence analysis around circular cylinders, hence, the most relevant studies on modeling of flow around circular piers are reviewed in the following, however, reports on morphodynamical simulations are barely considered here. The study of Olsen and Melaaen (1993) was one of the very first that dealt with three-dimensional flow simulation around a cylinder. Their model solved the Reynolds Averaged Navier-Stokes (RANS) equations on a structured, curvilinear grid, while turbulence was modeled with a k-ε turbulence closure. Using the results from flow calculations, the bed shear stress field was estimated together with an equilibrium sediment concentration, using an empirical formula. Based on those, local bed changes were calculated and a new iteration step started for flow calculations. Fairly good agreement between measured and calculated scour geometry was obtained, however, the calculated flow field was not verified against measurements. Richardson and Panchang (1998) also used a three-dimensional computational model (Flow-3D) to calculate the complex flow field, but in this case neither sediment transport nor scour formation was simulated. Detailed flow analyses were carried out for rigid flat bad case and for equilibrium scour geometry, too. They suggested the use of computational models contrary to empirical formulas to estimate scouring, however, the computational resources were stated as a limiting factor for such studies. Roulund et al. (2005) studied the flow field around a vertical circular pile exposed to a steady current. They used a RANS model with a k-ω turbulence closure, with no free-surface facility, to study the horseshoe vortex and lee-wake vortex flow processes around the pile. Furthermore, the influence of boundary-layer thickness, the Reynolds number and bed roughness on the horseshoe vortex was analyzed. The model validation was carried out against experimental data. Both steady and unsteady flow were simulated and they concluded that the timeaveraged bed shear stress obtained from the unsteady solution was found to agree better with measured data than the one obtained from the steady solution. This is important for the calculation of scour processes, where the computational time would be significantly larger for unsteady simulation. The flow model was coupled with a morphologic model to calculate scour around the cylinder and results showed that the simulation captured all the main features of the scour process. Numerical modeling of fluid flow has also been carried out in many other areas of water engineering. Chau and Jiang (2001) computed water quality in an estuary, while Ozmen-Cagata and Kocaman (2011) estimated the propagation of a dam break wave. Haun et al. (2011) solved the Navier-Stokes equations using a finite volume method to find the capacity of a spillway. In the last decade, due to the intensive growth of computer capacity and parallelized numerical modeling solutions, the development and application of large-eddy simulation (LES) became a convenient tool for the calculation of flow patterns around single and multiple cylinders, which result in very complex unsteady, threedimensional flows that are dominated by largescale coherent vortices (e.g. Palau-Salvador et al., 2008 and 2010; Stoesser et al., 2008 and2010). Tseng et al. (2000) carried out a detailed analysis to reveal the spatial flow structure around square and circular piers. Their results were compared with Dargahi's measurements, and good agreements were reported. They pointed out some slight differences between the strength of downflow, the horseshoe vortex as well as the wake vortex for the two obstacle geometries. Conclusions for the scour holing mechanism were also drawn, however, no morphodynamical calculations were performed. Catalano et al. (2003) performed LES to compute supercritical flow around a circular cylinder. To avoid the use of a high grid resolution at the walls a simple wall stress model was employed. They predicted the delayed boundary layer separation and drag coefficients which were consistent with measurements. Furthermore, their results were compared with steady and unsteady RANS simulations, and Catalano et al. (2003) provided evidence that RANS is generally less accurate than LES for such unsteady flows. LES together with laboratory-flume visualization were used to analyze coherent flow structures around a circular cylinder in a scour hole in the study of Kirkil et al. (2008). The bathymetry corresponded to the equilibrium scour case and was fixed in the simulations. A detailed investigation of the flow structure showed that the horseshoe vortex system comprised a primary necklace vortex along with several secondary vortices. They pointed out that the level of turbulent kinetic energy in the primary necklace vortex system is much larger than the one inside the approaching boundary layer. Their study also identified a mechanism for sediment entrainment in the regions of the two sides of the cylinder, induced by detached streaks of vorticity from the main necklace vortices.

EXPERIMENTS
The experiments that are used herein for validation purposes were carried out in the hydraulics laboratory of the Georgia Institute of Technology in a rectangular 24.4 m long, 1.1 m wide, and 0.5 m deep tilting flume. The water depth was controlled by an adjustable tailgate, which, in combination with the tilting mechanism, was used to produce uniform flow for a given flow rate. The test section, where the piers were located, was 2.75 m long and was filled with 0.18 m of sand. The parameters held constant for each case are the pier diameter, D = 0.0476 m, the flow rate, Q = 0.057 m 3 /s, and the uniform flow depth h = 0.1622 m. These conditions yielded the following dimensionless numbers, a Froude number of Fr = 0.26, a channel Reynolds number of Re = 53500, a friction Reynolds number of Re τ = 3900 and a pier Reynolds number of Re D = 15700. In total five experiments were conducted from which three were chosen for numerical model testing, as detailed flow measurements were available only for these cases. Before the experiments, with piers in place, were conducted, uniform flow condition was determined. The slope (S 0 ) of the flume for the experiments was set to S 0 = 0.0003639. Because the ultimate goal of this research was to investigate clear water scour around multiple piers, the critical velocity, which is defined as the flow velocity at which initiation of motion of sediment occurs, was estimated by Keulegan's equation (1938) to decide on the discharge (Q) which yielded Q = 0.057 m 3 /s. Discharge, slope and roughness conditions then yielded a uniform water depth, found by adjusting the tailgate. Flow uniformity for the flume without piers in place was then confirmed by water surface profile measurements along the flume centerline using a point gage with a precision of 0.0003 m. The uniform flow depth along the flume for each of the experiments was 0.162 m. For each experiment, pier 1 was located 17 m from the flume entrance. For case A, a 0.0476 m diameter PVC tube was used and was located on the centerline of the flume. For cases B and C two 0.0476 m diameter PVC tubes were used. The piers were spaced 0.191 m or 4D center to center, and were placed in the center of the flume. For the hydrodynamic results reported in this paper polyurethane was sprayed on the bed around the piers so that no scouring occurred. Velocity measurements were taken at selected verticals in the flow using Acoustic Doppler Velocimetry (ADV). For case A velocity measurements were taken along the centerline of the pier. For cases B and C a grid of 30 (x, y) locations with 8-10 points at fixed distances above the bed was taken.

Solver
The three-dimensional flow model solves the RANS equations discretized on a curvilinear grid system, built up by several structured blocks (Olsen, 2010) using the finite volume method. The discretization of the Navier-Stokes equations in the non-orthogonal grid follows the method given by Melaaen (1992). To calculate turbulence characteristics the model uses different zero-, oneand two-equation models, like Keefer's model (1971), the Spalart-Allmaras model (1992), the kω and the k-ε closure. For the latter one, the Kato-Launder modification (1993) can also be applied. In this study the well-known k-ε model was used to model the Reynolds stress terms. Using this assumption the Reynolds-averaged Navier-Stokes equations can be written as: where U = time-averaged velocity, t = time, ρ = fluid density, P = pressure, x j = Cartesian space co-ordinates, δ ij = Kronecker delta, k = turbulent kinetic energy and ν T = eddy viscosity coefficient. In the k-ε model ε represents the rate of turbulent energy dissipation. A transport equation is solved both for k and ε as a result of which the eddy viscosity coefficient can be evaluated as: The transport of k is modeled by the following differential equation: where P k defines the production of k, and this term is expressed as where S is the strain-rate, given as: Kato and Launder (1993) introduced a modification for this term in order to decrease the high turbulence viscosity in regions, where the normal k-ε model tends to produce too much turbulent energy, such as stagnation regions and in the small regions with strong acceleration and deceleration around the corners. With the modified production term they were able to produce much better results for transient simulations of vortex-shedding behind square cylinders. The proposal by Kato and Launder is to replace one of the strain-rates, S, in the turbulent production term with the vorticity, Ω. The Kato-Launder modified production then becomes: The transport of ε is modeled by the following differential equation: The constant values of the k-ε turbulence model are (Rodi, 1980 The pressure field is calculated using the SIMPLE algorithm (Patankar, 1980). In order to reach fast convergence, the free-surface facility was not used in this case, but a rigid lid approach was applied. As it was commented by Roulund et al. (2005) the effect of the water surface elevation difference between the upstream face and the sides of the pier on the velocity distribution in the vicinity of the pier is negligible for low Froude numbers. They showed that in case of a Froude number as high as 0.5 the water level difference can significantly influence the flow field, but not for lower values. In this study the Froude number was 0.26, thus the elimination of the free surface calculation is acceptable.

Boundary conditions
At the inlet, the streamwise velocities, U 1 were given corresponding to a log-law velocity profile. Transversal velocities, U 2 and vertical, U 3 velocity components were specified to be zero. For k and ε uniform distributions were defined with values typical for uniform flow conditions, using 0.0001 both for k and ε. At the outlet section, Neumann, i.e. zero-gradient boundary conditions were used for all variables. At the walls two different boundary conditions were used based on the dimensionless wall distance, defined as: where ν is the kinematic viscosity of the water. According to Schlichting's experiments (1979) the formula for rough walls is valid in the range where y + is between 30 and 3000. For lower values of y + the smooth wall-laws is to be used. For the bottom, where the bed material with d 50 = 0.0011 m was placed (where d 50 is the grain diameter for which 50% of the sediment sample is finer by weight) we applied the rough wall-laws: Here, U = velocity parallel to boundary layer, u * = friction velocity, κ = von Karman constant (κ =0.41), y = distance between wall and the investigated point and k s = roughness height, defined as 3d 90 according to van Rijn (1984), where d 90 is the grain diameter for which 90% of the sediment sample is finer by weight. For the smooth side walls and the cylinders, the following function was used: (13) where E is an empirical parameter equal to 9.0. For k and ε the formulas given by Rodi (1980) were defined, assuming that the turbulent production is equal to the dissipation of k near the wall. The formula for ε is given as: To find k at the wall, the source of the transport equation for k is integrated over the bed cell including the formula for ε. See the details in the note of Olsen (2000). Zero gradient boundary conditions were used at the water surface for the horizontal velocities and ε. The turbulent kinetic energy and the vertical water velocity were set to zero at the water surface. The SIMPLE method does not need boundary conditions for the pressure.

Nested grid algorithm
The nested grid approach can be used where the domain has important local flow features that feature different scales. The method consists of a fine grid placed inside a coarse grid. Typically, the coarse grid covers the river while the fine grid resolves the smaller flow features around a structure. The main numerical algorithms needed for the system are: 1. How the parameters in the coarse grid affect the values in the fine grid.
2. How the values in the fine grid affect the coarse grid flow field.
The first question is about interpolation of values from the coarse grid to the boundary of the fine grid. The fine grid requires some boundary conditions. The approach presented here is to look at each outer surface, s, of the cells in the fine grid. A search algorithm is first used to find the four cells in the coarse grid that are closest to s, but on the outside of the fine grid domain. These four cells were called c i , where i=1.. 4. The interpolation of the value of a parameter f on the surface s is then given from the following formula, where u c,i is the velocity vector in the coarse cell i: The vector pointing from the coarse grid to the center of surface s is denoted n s . The dot product is used to give more weight to the coarse cells upstream of the surface. If the dot product is negative from all the four coarse cells, then the surface has an outflow. Boundary conditions for a general variable are then not needed. This formula gives good results without any need for adjusting parameters. Question 2 is not so straightforward. If the constructions in the river are very small, then this effect may not be very important. However, for bridge piers with large diameters, the effect can be significant. Several different algorithms to interpolate the values from the fine grid onto the coarse background grid that is located underneath the fine grid are tested. However, this treatment is prone to testability problems and nonconvergence. An alternative is to ignore the coarse grid results of the background grid and only adjust the coarse grid computation on the outside border of the fine grid using a source term. The equation for the coarse cell, indicated by subscript p, closest to the surface area which bordered the fine grid cell, denoted with subscript f, is modified as: The * denotes the original values. The j parameter in the formula is the number of fine grid cells for each coarse cell. The parameter r is a relaxation factor set to 0.01. The method is in principle independent of the alignment of the fine and the coarse grid (Olsen and Tjomsland, 1998). However, the method works best in areas where the effect of the fine grid flow field is less significant on the coarse grid parameters.

Numerical grids
The computational grids consist of two blocks in all cases. The channel is discretized using a coarse, orthogonal grid with a length of 3 meters with a horizontal resolution of 0.05 m, and consisted of 12 layers in the vertical direction, which yield 60 x 20 x 12 cells. A much finer second block, the nested grid is implemented in the coarse grid in the surroundings of the cylinder(s). In the nested grid, exploiting the flexibility of the curvilinear description, the grid resolution is substantially increased, particularly at the cylinder wall(s), to obtain a more accurate estimation of the flow field, especially to enable capturing accurately the flow separation (see

RESULTS
Two different data sets are available from the laboratory experiments. For the single pier case four profiles were measured along the flume centerline. One upstream of the pier at a distance of three-diameter and three downstream of the pier at one, two and three-diameter distances from the pier. A comparative analysis of the streamwise velocity and the turbulent kinetic energy (TKE) profiles is carried out. More detailed flow measurements were done for the multiple pier cases, so that three longitudinal 2D planes of flow velocities and TKE along the centerline of pier 1, the centerline of the flume and the centerline of pier 2 are available for comparison.

Single pier case
As a first step of the numerical investigations a sensitivity analysis was carried out looking at the influence of horizontal and vertical grid resolution, nested grid interpolation parameters, time step, wall laws etc. When using smooth wall laws the nested grid resolution close to the piers apparently influences the flow separation. In order to describe the flow separation accurately three grid resolutions are tested with an average grid size of 0.005 m, 0.0025 m and 0.0015 m. For the first two cases the flow separation is poorly reproduced, but defining a grid resolution of 0.0015 m, moreover, refining the grid at the pier, a realistic flow pattern is achieved. The test suggests the use of the finest mesh, however, the computational demand obviously increases. The computational time on an IBM p575+ computer with 16 cores and a processor speed of 1.9 GHz is approximately 40 minutes. Fig. 5 shows randomly distributed velocity vectors close to the water surface for the three grid resolutions. The use of 12 layers vertically yielded an optimal resolution, since there is no significant change in the flow pattern when employing a finer resolution of 20 layers, as introduced in the following. Fig. 6 shows U 1 profiles for the four measurement locations indicating the calculated values with and without the Kato and Launder modification. For the former case the model results with 20 vertical layers are also plotted in order to assess the effect of vertical grid resolution. The approach flow boundary layer is reproduced satisfactorily. The measured velocity profiles downstream of the pier suggest complex flow. As a result of the horseshoe vortex system, higher streamwise velocities characterize the near bed region, whereas the negative values close to the water surface indicate that the extent of the recirculation zone behind the pier is about 2D. The results from the numerical simulation in downstream of the cylinder show somewhat smaller velocities close to the bottom than the experiments, i.e. the recirculation zone is slightly overpredicted here. Near the water surface, the agreement is quite good. The disagreement could be a result of the isotropic turbulence assumption, but on the whole, the agreement is satisfying. The velocity profiles from the 20 vertical layers case have almost the same shape, however, a slight difference at the last point can be observed. It can be stated that the higher vertical resolution has no major influence on the flow. As to the turbulence profiles (Fig. 7), the approaching flow is characterized by a reasonably low and uniform turbulence level. The numerical model predictions are in pretty good agreement with the measurements. Downstream of the cylinder the level of turbulent kinetic energy increases by approximately one order of magnitude. Surprisingly, in 1D distance downstream of the pier the measured turbulence profile shows more uniformity than the one from RANS simulation. As expected from previous studies (e.g. Kato and Launder, 1993) the k-ε   With the Kato and Launder modification the turbulence is better captured. The correction influences only slightly the velocity distribution (see Fig. 6), but at the same time the increase in turbulent kinetic energy towards the water surface is well predicted for the three profiles downstream of the pier. Particularly the two last profiles agree very well with the measured ones.
It should be noted that the underestimation of turbulent energy close to the flume bottom for the second profile will result in lower bed shear stress and would lead to inaccuracies in the morphological modeling. Similar to the velocity profiles, the turbulence profiles also suggest that higher vertical resolution has no significant effect on the flow results. The three-dimensional nature of the flow structure can be assessed by looking at the vertical velocity components in the vicinity of the pier.  Dargahi (1989). Similarly, smaller vortices are generated along the pier wall and also behind the pier. Large-scale vortices are found in the upper region downstream of the pier, the size of which corresponding to the pier diameter.

Double-pier cases
In order to validate the numerical model, distributions of U 1 , U 3 and TKE distributions in longitudinal planes are compared. To enable easy   comparisons model results are only shown below 0.10 m. Though the magnitude of the streamwise velocity is slightly overpredicted along the centerline for case B (Fig. 11), the spatial variation of the velocity shows a very similar pattern to the one in the experiment. Along the centerline of pier 1 the size and shape of the recirculation zone behind the cylinder is well predicted, moreover, increased streamwise velocity near the bed further downstream is observed, matching well with the experiment. Also, the boundary layer flow approaching pier 2 agrees well with the measurements. As to the vertical velocity components (Fig. 12), the calculated strong upwelling behind pier 1 is predicted well and has a similar magnitude as in the measurements. Along the centerline of the  channel, i.e. between the two piers, a weak downwelling is computed, suggesting secondary flow there. The strength of the secondary current is slightly over-predicted by the RANS model. The model predicts strong downwelling in front of pier 2, initiating the horseshoe vortex system (e.g. Roulund et al., 2005). A uniform distribution of TKE characterizes the flow along the centerline and the line of pier 2 throughout the water depth except close to the bed, where higher values of TKE are observed due to the high shear stress there (Fig. 13). There is a  strongly turbulent patch at the upstream toe of pier 2 indicating the influence of downwelling and the initiation of the horseshoe vortex there and the reason for the start of the scour hole development around such obstacles. This is an important feature to be predicted by numerical models for sediment transport and bed morphology development. Similarly to the single pier case, the turbulent energy suddenly increases behind pier 1. Despite the time-averaged description of the flow the effect of the coherent vortices is reflected in the turbulence distribution, and this is in accordance with the laboratory experiments. As the plots of U 1 show, the length of the recirculation zone is approximately 2D, in which less turbulence is expected. Indeed, an increase of TKE is found after 2D and due to the upwelling the maxima are found close to water surface. This is in good agreement with the experiments.
The three-dimensional flow structure can be visualized through pathlines starting from distinct points upstream of the piers (Fig. 14). Particles released close to the bed and in the centerline of the piers will enter the horseshoe vortex system having a downwelling character in front of and around the cylinders. Behind the pier they travel in a spiraling motion towards the water surface.
There are no significant differences compared to the single pier case, i.e. the flow pattern is the result of the individual piers. Hence, the interaction between the piers is not significant. This behavior, however, may not be the same for other pier layouts, especially if the obstacles are located closer to each other. Similar to case B the streamwise (Fig. 15) and the vertical velocities (Fig. 16) are well reproduced in case C. The streamwise velocities in the centerplane are somewhat smaller than those in case B, which is due to less flow contraction. Although the distance between the piers is the same as in case B, the distance projected onto the plane perpendicular to the main flow direction is larger. The measured vertical velocities close to the water surface along the centerline of the channel indicate permanent downward motion suggesting secondary flow there. The downwelling is reproduced by the model, however, slightly underestimated. Contrarily, the upward motion of the flow behind pier 1 is overestimated to some extent.   Calculated TKE fields (Fig. 17) agree reasonably well with the ones from the experiments, the highest values are under-predicted, though. The interaction between the piers is even weaker than that in case B, and unstable vortices dominate the flow behind pier 1, but are less influenced by the flow through the piers resulting in higher deviation of the velocities, i.e. turbulent energy. The numerical model results in more intensive flow contraction yielding (as seen in Fig. 15) less turbulent flow conditions here.
A key input of morphological modeling is the distribution of the bed shear stress as calculated by the hydrodynamic model. Using the k-ε model the bed shear stress value (τ) is commonly estimated based on the assumption that the production and dissipation of turbulence is in equilibrium close to the wall. Hence, the bed shear stress can be computed directly from the k value near the bed: where C μ = constant for the k-ε model (introduced above), and ρ = density of water. In order to see if the calculated bed shear stress field is acceptable for morphodynamical investigations, the calculated and measured values are compared. However, direct bed shear stress measurements were not available from the laboratory experiments. Therefore, two methods are used to estimate the bed shear stress: i) based on the measured turbulent kinetic energy closest to the bottom (measured in 0.005 m distance from bed) for which Eq. (18) is used; ii) the logarithmic velocity profile (Eq. 10) to be fitted to the three lowest measurement points using the leastsquares method, from which u * is calculated. From u * the bed shear stress is computed with: Longitudinal profiles of bed shear stress distributions are plotted along the centerline of pier 1, along the flume centerline and along the centerline of pier 2, ( 60% higher than in case B), where the calculations result in under-prediction due to the time averaged description. The same underprediction was introduced for the single pier case. The accurate estimation of bed shear stress in a considerably unsteady flow pattern is a limitation of RANS models. Nevertheless, due to the relatively low computational demand of RANS models, they are a suitable tool for practical river engineering problems, where local flow and morphological features are to be studied. As a reference, the computational time (IBM p575+ computer with 16 cores and a processor speed of 1.9 GHz) for case B and case C were 90 minutes and 140 minutes, respectively.

CONCLUSIONS
A three-dimensional numerical model was developed, which solves the Reynolds Averaged Navier-Stokes equations with a k-ε turbulence closure on a curvilinear nested grid system. For model validation a laboratory experiment was chosen, in which the flow and turbulence field were studied for a single-and for two double circular cylinder cases. The following conclusions could be drawn: 3. The calculated three-dimensional flow structure showing the strong downflow in front of the pier, the spiral upwelling vortices behind the pier together with the relatively high TKE patches indicates the existence of the horseshoe vortex system. It proves that the RANS description might be able to reproduce the complex flow field in such a case, however, very fine grid resolution was needed in the vicinity of the pier in order to capture correctly the location of flow separation.
4. The extracted longitudinal distributions of streamwise and vertical velocity components from the numerical model predictions for the double pier cases are in accord with measured distributions. The characteristic features, i.e. a turbulent boundary layer of the approaching flow, flow contraction between the piers, flow separation and strong upwelling flow behind the piers, are predicted fairly well suggesting that the model is capable of simulating such complex flows.
5. Two methods were applied to calculate bed shear stress values from ADV measurements, firstly from measured TKE values close to the bed and secondly using the logarithmic law of the wall, in which measured velocities are fitted to the log law which then yields the bed shear velocity. In the undisturbed boundary layer both approaches give similar bed shear stress values, and the numerical model predictions are in good agreement with those. Behind the cylinder, the second method is not feasible as the logarithmic distribution is not valid in this zone. For case B the calculated bed shear stress distribution behind pier 1 shows almost perfect agreement with the values estimated from ADV measurements. In case C the numerical simulation under predicts the bed shear stress in the wake zone, when compared with the TKE obtained stresses. The same behavior is true for the single pier case, indicating that for case C, there is little interaction between the two piers. High values of bed shear stress are found in front of all the piers indicating potential for scouring here.
6. Using the nested grid algorithm the threedimensional RANS model is capable of predicting the complex flow and turbulence features in the vicinity of structures obstructing a flow. The computational demand of such a simulation is significantly lower than discretizing the whole domain with a high resolution grid. This is important for the numerical modeling of local scouring, which becomes feasible using such a tool.
Though unsteady features such as vortex shedding cannot be captured by a RANS approach and hence certain inaccuracies are to be expected, the introduced method is applicable to real river engineering studies in which local complex flow features are investigated.