1Introduction

Flows past stationary structures in channels, rivers, harbors, and coastal areas are of considerable importance because of sediment transport and scouring which results in failure of the structures. Understanding of the behavior of these flows and bottom sediment transport requires a series of experimental investigations. However, physical experiments demand an enormous time and include a similarity limitation between the prototype and the model (Ettema et al., 2004). Therefore an alternative method is desirable. Thus many recent studies have been devoted to the development of numerical modeling.

Simple attempts to model the nature of turbulence in the wake of flow past stationary objects by 2-D numerical simulations often use cross section-averaged or depth-averaged flows(Bosch and Rodi, 1998; Park et al., 1998; Chen et al., 2003; Honzejk and Frana, 2008). Numerical modeling allows us to elucidate the comprehensive physical mechanism of the wake. Olsen and Melaaen(1993) performed a 3-D numerical simulation to estimate the turbulent flow field around a cylinder structure. Ishino et al. (1993) examined wake characteristics behind a cylinder as the Reynolds number increased from small to large values. Richardson and Panchang(1998) produced turbulent motion near the structures and compared model results with experimental data. To investigate wakes and scouring problems, some studies applied a computational fluid dynamics code with a built-in turbulence model, FLOW-3D(Yang and Choi, 2002; Park et al., 2009).

Wakes behind a structure depend on the drag coefficient of the structure and the flow characteristics. Accurate predictions of wakes are critical to the safety of structures in channels, rivers, harbors and coastal areas. The objective of this study is to calculate wakes behind a square cylinder using two kinds of different turbulence models based on the eddy viscosity concept such as the zero- and the one-equation model. In addition, we evaluated each model’s performance by comparing its results with analytical solutions using three skill assessments: the correlation coefficient(r) for the similarity of the wake shape, the error of maximum velocity difference(EMVD) for the accuracy of wake velocity and the ratio of drag coefficient(RDC) for the pressure distribution around the structure. It was also investigated the applicable range of empirical constants and suitable values for each turbulence model.

2Analytical Solutions and Turbulence Models

2.1Analytical Solutions

Schlichting(1930) investigated wakes behind a single body based on the Prandtl’s mixing length hypothesis and suggested an equation for a wake. The solution of the equation can be written as follows:

where U is the free stream velocity; u_D is the velocity difference of a wake(U - u); u is the flowing velocity; l is the mixing length; x is the flowing distance; y is the transverse distance; C_D is the drag coefficient of a resistance body; D is a representative distance of the body such as the diameter; b is half of the width of a wake; and β is the rate of l and b, which can be determined by measurement values.

Görtler(1942) and Reichardt(1951) introduced the eddy viscosity ν_t that is a function of the non-zero mean velocity gradient, and assumed that ν_t is a constant ν₀ . They provided an equation for a wake using ν_t and its solution can be given by:

These theoretical solutions, Eq.(1b) and Eq.(2b), are comparable and have been compared with Schlichting’s measurements (Schlichting, 1979). There is excellent agreement between the theoretical solutions and measurements, as shown in Fig. 1. Therefore the solution especially Eq.(2b) was used to evaluate wakes behind a square cylinder in this study.

2.2Turbulence Models

Turbulence models are based on the turbulent-viscosity hypothesis introduced by Boussinesq in 1877, which states that the Reynolds stress is proportional to the mean strain rate. As a result, the eddy viscosity can be specified as a function of the characteristic velocity and the characteristic length. In the mixing length or zero-equation model, the characteristic velocity is determined by the mean velocity gradient; thus, the eddy viscosity is obtained as:

where ∂ u / ∂ y is the mean velocity gradient. It is wellknown that the zero-equation model works reasonably well in case of simple turbulent shear flows. However the applicability of the model is limited when turbulence transport becomes small in a flow in which the production P and dissipation ε of turbulent kinetic energy are approximately in balance(P / ε ≈ 1).

The one-equation model employs turbulent kinetic energy k to calculate the characteristic velocity. In the FLOW-3D model(Flow Science, 1999), the k is governed by:

where t is time; V_F is the fractional volume; u, v, w are the velocity components in the x-, y- and z-direction, respectively; A_x , A_y , A_z are the fractional area open to flow in the x-, y- and z-direction, respectively; G is a buoyant production; and D_k is a diffusion term. Kolmogorov(1942) and Prandtl(1945) suggested that the characteristic velocity is proportional to k and can be described as: where L is the empirical mixing length scale. Wilcox(1993) noted that the one-equation model has an advantage in accuracy over the zero-equation model based on a comparison of model predictions with measurements.

3Numerical Models and Simulation Conditions

3.1Numerical Models

The FLOW-3D model based on the marker/cell and the solution algorithm volume of fluid(SOLA-VOF) methodologies was used in this study. The two methodologies were developed by Hirt and Nichols(1981) and applied to a wide range of complicated problems(Flow Science, 1999). To calculate unsteady flows near the areas of the fluid/obstacle interface, we used the FLOW-3D using the fractional area in which the volume-obstacle representation method incorporated into equations of motion discretized by the finite volume method. The model uses a non-uniform staggered grid so that velocity and stress are defined at the cell face, and the scalars are located at the cell center. The velocity and pressures were computed using the continuity and momentum equations, and the water elevations were solved using the VOF method. To calculate turbulence flow around a square cylinder, we utilized the built-in turbulence modules in the FLOW-3D. A detailed description of the model is given the FLOW-3D User’s Manual(Flow Science, 1999)

3.2Drag Coefficients and Correlation Coefficients

When a flow moves past a blunt object, the interaction between the object and the flow can be described as the drag force obtained by the averaged pressure forces, neglecting the skin friction of the object. The averaged pressure forces are calculated by integrating the x-, y- and z-component of the pressure force. The drag force can be represented as a dimensionless drag coefficient, and the drag coefficient in the x-direction is defined as:

where F_x is the drag in the x-direction; A_x′ is the projected area in the x-direction; and U_x is the upstream velocity in the x-direction. The drag coefficients in the y- and z - directions are similar to the coefficient in the x -direction. The drag force at each cell in the FLOW-3D was obtained as a product of the pressure and projected area in the x-, y-, and z - direction such as Eq.7(a), and can be separated into the x-, y-, and z - component of the drag force. The total drag force of the object can be written as Eq.(7b): where p is the pressure at a cell; dA′ is the projected area in the x-, y-, and z -direction at a cell; F is the drag in a cell; and i, j, k are unit vectors in the x -, y-, and z -direction, respectively(Flow Science, 1999).

The standard correlation coefficient r(Kang et al., 1993) shows the statistical relationships of two variables and is bounded from -1 to 1. In this study, the r approaches 1 as the agreement between analytical and simulated velocities would be increased. A good agreement between two velocities means that the model successfully represents the wake shape and width.

3.3Simulation Conditions

For adequate calculation of wakes behind a square cylinder, we used a numerical basin 280 m long, 60 m wide, and 5 m deep (Fig. 2). A non-uniform grid containing 190×44 horizontal cells and 11 vertical layers was used to resolve the region of high flow gradients shown in Fig. 2(b). The range of grid sizes is from ∆x × ∆y × ∆zmin = (0.333m × 0.333m × 0.5m) to ∆x × ∆y × ∆zmax = (2.641m × 2.508m × 1.0m). A constant flow propagates toward a square cylinder that is 30 m from the boundary inlet. The size of a free-standing obstacle is (1 m × 1 m × 5 m) = (length × width × height) and the drag coefficient of the square cylinder is 2.2 at a Reynolds number Re ~ 1.0 × 10⁵(Blevins,1984 ).

At the boundary inlet, a constant velocity, V = (u , 0, 0) = 0.1 m/s, was specified. The two side boundaries were specified with a symmetry condition. The continuative boundary condition was applied at the flow outlet. Meanwhile the bottom boundary conditions were as follows: the smooth bottom, the no-slip condition, and the normal velocity equals zero at the bottom.

The FLOW-3D with inclusion of the zero- and the one-equation model was used to simulate wakes behind the square cylinder, in which these turbulence models included an empirical constant such as the mixing length. In this study, several wake simulations(Table 1) were conducted with values of the mixing length from 0.3b_1/2 to 0.8b_1/2 increasing by 0.1b_1/2 which was employed from the similar wake experiments of Rodi(1972) and Schlichting(1930). We used a Reynolds number of 1.0 × 10⁵ that was adopted from the Blevins(1984)’s results. In addition to remove the transient effects of the initial condition, the simulation was performed for over 5000 seconds with varying time-steps controlled by the program. After the flow reaches a steady state condition and the simulation results were analyzed.

4Model Results and Discussion

4.1General Horizontal and Vertical Distribution of Wakes

Using the one-equation model, we calculated horizontal velocity fields and longitudinal(x-direction) velocity along the centerline at the surface, middle, and bottom layers(Fig. 3). As can be seen in Fig. 3(a), wakes behind the square cylinder were similar to the long tail of a shooting star in the sky and the width of wake increased as the flow moved downstream. The magnitude of the horizontal velocity at the surface layer was larger than that of the middle and bottom layers, leading to a wide wake at the surface layer. From the longitudinal velocity along the centerline, which was non-dimensionalized using the constant velocity of the basin, we note that a reverse flow occurred near the square cylinder due to the low pressure. After passing the square cylinder, the flow recovered rapidly. Where the value of x/ (C_DD) was greater than 10, the velocity gradient slowly varied and became stabilized after the value of x/ (C_DD) was greater than around 50. These results are quite close to those obtained by Schlichting(1930). Horizontal velocity at the surface layer was generally larger than at the other layers. Resulting from strong vertical mixing and interaction of the vortex, the differences in horizontal velocity at the three layers are insignificant just behind the structure. On the other hand, the results of the zero-equation model were similar with those of the one-equation model, so not presented in here.

4.2The Zero-Equation Model

In the zero-equation model or mixing length model, the eddy viscosity becomes a product of the mixing length l and the velocity scale l∂ u / ∂ y . The model is arguably the simplest turbulence model; its major drawback is that the mixing length must be empirically specified. For a complex flow, the appropriate specification is very difficult and requires considerable guesswork. Consequently, we have little confidence in the accuracy of the results. Various values of the mixing length were employed to represent the wake velocity and width behind the square cylinder(Fig. 4). For a small value of the mixing length l = 0.3 b_1/2, model results underestimated the wake width but its width increased as the mixing length increased. When the mixing length was close to 0.5 b_1/2 the wake width from model simulation and the analytical solution were in good agreement.

Using the correlation coefficient(r), the error of maximum velocity difference(EMVD), and the ratio of drag coefficient (RDC), it was evaluated the performance of the mixing length model for the simulation of wakes behind the square cylinder (Table 2) where all r, EMVD, and RDC values were space averaging values. The values of r between the mixing length model predictions and the analytical solutions ranged from 0.954 to 0.997, thus indicating significant correlations as shown in Table 2. The minimum EMVD of 0.426 occurred when l was equal to 0.4b_1/2 . The values of RDC were 0.965 and 1.198 for l = 0.3b_1/2 and l = 0.4b_1/2, respectively; thus, a perfect RDC value of 1 could be obtained between l = 0.3b_1/2 and l = 0.4b_1/2 . From the test cases, the absolute mean error(AME) integrating r, EMVD and RDC ranged from 20.3 % to 56.3 %, which was relatively high compared to other model results. In particular, the model overestimated the velocity difference, whose EMVD ranged from 42.6 % to 66.0 %. Schlichting(1930) suggested a value of l = 0.41b_1/2 ≈ 0.4b_1/2 by a theoretical study. However the model simulation with this value created a high EMVD of 42.6 %, resulting from limitations inherent in the model which underestimate the eddy viscosity near the velocity gradient of zero behind the square cylinder centerline. In the case of plane jet flow, Rodi(1972) used the l value of 0.254b_1/2 but this value is not suitable for the simulation of wakes behind the square cylinder based on the model results at l = 0.3b_1/2 . As a consequence the mixing length model is not applicable to calculate wakes behind the square cylinder.

4.3The One-Equation Model

The one-equation model uses the k instead of the mean velocity gradient to calculate the characteristic velocity because the velocity becomes zero where the mean velocity gradient is zero in the zero-equation model. There are some cases in which the mean velocity gradient is zero but the eddy viscosity is non-zero. Thus it is better to base the velocity scale on k and the one-equation model has an advantage in accuracy over the zero-equation model(Wilcox, 1993). There remains a problem that the mixing length scale L must be specified even though the velocity scale is obtained by solving the k-transport equation. The estimation of L for the simulation of wakes behind the square cylinder is a challenge because of the uncertainty of L. To acquire the appropriate value of L, we simulated wakes using the one-equation model with a variety of L values as used the same range of the mixing length in the tests of zero-equation model(Table 1).

The same as the zero-equation model results, the width of wake increased as the L increased(Fig. 5). The range 0.952 to 0.997 of r shows significant correlation between the model predictions and analytic solutions. In particular, the values of r were near 1 (0.993 to 0.997) when the L range was 0.5b_1/2 to 0.7b_1/2 and within these ranges the width of the wake barely changed. The EMVD indicated that the model underestimates the velocity difference, which is the opposite of the zero-equation model results. The range of EMVD was -5.9 % to -34.2 %, which is relatively low compared to the zero-equation model, and the minimum EMVD occurred when the L was equal to 0.5b_1/2. The value of RDC ranged from 0.7 to 1.245 with an almost perfect fit(≈ 1.002) at L = 0.6b_1/2and the RDC value increased as the L increased. The value of AME ranged from 3.4 % to 19.9 % and, in particular, when the L was between 0.5b_1/2 and 0.6b_1/2, the value of AME was less than 5.6 % as shown in Table 3. On the basis of comparisons between the zero- and the one-equation model, the eddy viscosity calculated using turbulent kinetic energy better represented wakes behind the square cylinder than the eddy viscosity obtained using the mean velocity gradient. An L value of 0.6b_1/2 in the one-equation model seemed to be the best value to simulate wakes behind the square cylinder.

5Conclusion

It was simulated wakes behind a square cylinder using two different turbulence models based on the eddy viscosity concept, such as the zero- and the one-equation model and investigated the applicable range of the major empirical constant in each model. Based on the results, the following conclusions are drawn: