• For Contributors +
• Journal Search +
Journal Search Engine
ISSN : 1229-3431(Print)
ISSN : 2287-3341(Online)
Journal of the Korean Society of Marine Environment and Safety Vol.25 No.3 pp.344-353
DOI : https://doi.org/10.7837/kosomes.2019.25.3.344

# Analysis of Empirical Constant of Eddy Viscosity by k-ε and RNG k-ε Turbulence Model in Wake Simulation

Il Heum Park*, Young Jun Cho**, Jong Sup Lee***
*School of Marine Technology, Chonnam National University, Jeonnam 59626, Korea
**Ocean and Port Research, Jeonnam 59643, Korea
***Department of Civil Engineering, Pukyong National University, Busan 48513, Korea
Corresponding Author : parkih@chonnam.ac.kr, 061-659-7152
May 13, 2019 May 24, 2019 May 28, 2019

## Abstract

The wakes behind a square cylinder were simulated using two-equation turbulence models, k-& and RNG k-& models. For comparisons between the model predictions and analytical solutions, we employed three skill assessments:, the correlation coefficient for the similarity of the wake shape, the error of maximum velocity difference (EMVD) of the accuracy of wake velocity, and the ratio of drag coefficient (RDC) for the flow patterns as in the authors’ previous study. On the basis of the calculated results, we discussed the feasibility of each model for wake simulation and suggested a suitable value for an eddy viscosity related constant in each turbulence model. The k-& model underestimated the drag coefficient by over 40 %, and its performance was worse than that in the previous study with one-equation and mixing length models, resulting from the empirical constants in the &-equation. In the RNG k-& model experiments, when an eddy viscosity related constant was six times higher than the suggested value, the model results were yielded good predictions compared with the analytical solutions. Then, the values of EMVD and RDC were 3.8 % and 3.2 %, respectively. The results of the turbulence model simulations indicated that the RNG k-& model results successfully represented wakes behind the square cylinder, and the mean error for all skill assessments was less than 4 %.

## 초록

Chonnam National University

## 1. Introduction

An understanding of the behavior of turbulent flow and bottom sediment transport behind structures 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. Many studies have been devoted to the development of numerical modeling in the past decades (Bosch and Rodi, 1998;Liou et al., 2002;Chen et al., 2003;Honzejk and Frana, 2008;Park et al., 1998;Park et al., 2009;Park et al., 2014).

Simple attempts have been made to model the nature of turbulence in the wake of the flow past the stationary objects engaged in 1-D or 2-D numerical simulations using the cross section averaged or depth averaged flows (Bosch and Rodi, 1998;Park et al., 1998;Honzejk and Frana, 2008). With computer hardware and numerical scheme developments, numerical modeling allows us to elucidate the comprehensive physical mechanism of the wake (Park et al., 2014). In previous studies, Olsen and Melaaen (1993) performed a 3-D numerical simulation with the k-ε model to estimate the turbulent flow field around the cylinder structure, and Ishino et al. (1993) examined wake characteristics behind the cylinder as the Reynolds number increased from small to large values. The turbulent motions near the structures were produced and compared with the experimental data from Richardson and Panchang (1998) who utilized the RNG (Renormalized Group) k-ε model. To analyze the flow around the square and circular structures, the Large Eddy Simulation (LES) turbulence model was adopted by Tseng et al. (2000), Liou et al. (2002) and Chen et al. (2003). Yang and Choi (2002) simulated the turbulent wake behind the bridge piers with the built-in RNG k-ε and LES models of the computational fluid dynamic commercial code (FLOW-3D). Park et al. (2014) recently simulated wakes behind a square cylinder using the zero- and one-equation turbulence models with FLOW-3D and investigated the applicable range of the major empirical constant in each model.

The wakes behind a structure depend on the drag coefficient and the flow characteristics. Accurate predictions of wakes are critical to the safety of structures in channel, rivers, harbors and coastal areas (Park et al., 2014). The objective of this study is to calculate the wakes behind a square cylinder using two types of two-equation turbulence models based on the eddy viscosity concept, the k-ε and RNG k-ε models. In addition, we evaluate each model performance by comparing its results with analytical solutions using three skill assessments: the correlation coefficient (γ) 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 (see Park et al., 2014). The study also investigates the applicable range of empirical constants and suitable values for each turbulence model in the same manner as Park et al. (2014).

## 2. Analytical Solutions and Turbulence Theories

### 2.1 Analytical Solutions

Schlichting (1930) investigated wakes behind a single body based on Prandtl’s mixing length hypothesis and suggested an analytical solution. 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 ν0. They provided an equation for a wake using νt and its solution can be given by

$u D U = 1 4 π U C D D ν 0 x C D D − 1 / 2 exp − y 4 U ν 0 x$
(1)

where uD is the velocity difference of a wake, U is the free stream velocity, CD is the drag coefficient of a resistance body, D is the representative distance of the body, x is the flowing distance, y is the transverse distance and νt = ν0 = 0.0222CDD U . This theoretical solution Eq. (1) was in excellent agreement with Schlichting’s measurements (1979) and was used to evaluate, therefore, the wakes behind a square cylinder in this study (see Fig. 1 of Park et al., 2014).

### 2.2 Turbulence Theories

For the broadest range of applicability of the turbulent model, the characteristic length scale L and the velocity scale need to be specified through theoretical approaches. In the one-equation model, the characteristic velocity scale was expressed as the function of the turbulent kinetic energy k as $ν t = L 2 k / 3$. The turbulence energy dissipation ε scaled as Ck3/2/L where C was a constant (Chou, 1945). Thus k and ε could be combined to eliminate the empirical mixing length and to produce an eddy viscosity equation Eq. (2). The k-ε model, one of the two-equation models, solved two transport equations for k and ε and, after that, yielded the eddy viscosity.

$ν t = C μ k 2 ∈$
(2)

where Cμ is an empirical constant and the typical Cμ is Cμ0 (= 0.09) (Launder and Spalding, 1974).

Using a field-theoretical approach and statistical methods, the RNG k-ε model was used to calculate model constants. The basic idea was that the small turbulence scales could be removed through space-time Fourier decomposition of the velocity field and repeated substitution in the Navier-Stokes equations. Also, the model constants were evaluated from the large scale field and the modified viscosity induced by the elimination of the small scales (Yakhot and Orszag, 1986;Yakhot et al., 1992). The model introduced the modified eddy viscosity νeddy in lieu of an eddy viscosity in the k-ε model (Yakhot et al., 1992).

$ν e d d y = ν 1 + C μ R ν k ∈ 2$
(3)

where ν is the kinematic viscosity and CμR is an empirical constant and its typical value is known as CμR0 (= 0.085) (Yakhot and Smith, 1992;Yakhot et al., 1992).

## 3. Numerical Models and Simulation Conditions

### 3.1 Numerical Models

FLOW-3D models (Flow Science, 1999) based on the Marker and Cell (MAC) and the Solution Algorithm Volume Of Fluid (SOLA-VOF) methodologies were used in this study. The velocity and pressures were computed using the continuity and momentum equations, and the water elevations were solved using the VOF (Volume of Fluid) method (Flow Science, 1999). To calculate the turbulence flow around a square cylinder, we utilized built-in turbulence modules such as the k-ε and RNG k-ε models in FLOW-3D. A detailed description of these models is given in the FLOW-3D User’s Manual (Flow Science, 1999).

The principal goal of any turbulence model is to provide a mechanism for estimating the influence of turbulent fluctuations on mean flow quantities. This influence is usually expressed by the additional diffusion terms in the equations for mean mass and momentum. In the one-equation or k-ε models, k is governed by Eq. (4).

$∂ k ∂ t + 1 V F u A x ∂ k ∂ x + υ A y ∂ k ∂ y + w A z ∂ k ∂ z = P + D k − ∈$
(4)

where P is a shear production and Dk is a diffusion term. The transport equation for the ε is

$∂ ∈ ∂ t + 1 V F u A x ∂ ∈ ∂ x + υ A y ∂ ∈ ∂ y + w A z ∂ ∈ ∂ z = ∈ C ∈ 1 k P + D ∈ − C ∈ 2 ∈ 2 k$
(5)

where C1 and C2 are the user defined constants and D is a diffusion term. In the RNG k-ε model (Yakhot and Smith, 1992;Yakhot et al., 1992), the additional term R on the righthand side of Eq. (5) is an extra strain term given by

$R = C μ R η 3 1 − η / η 0 1 + β R η 3 ∈ 2 k$
(6)

Where , η0 = 4.38 and βR = 0.012.

### 3.2 Drag Coefficients

As described in the previous study by Park et al. (2014), the drag coefficient in the x-direction is defined as

$C D x = F x 0.5 ρ A x ′ U x 2$
(7)

where Fx is the drag force in the x-direction, Ax′ is the projected area in the x-direction and Ux is the upstream velocity in the x-direction. Park et al. (2009) and Park et al. (2014) calculated the drag coefficient of an object in the numerical basin using FLOW-3D, and the same approaches are applied to this study.

### 3.3 Simulation Conditions

For the adequate calculation of wakes behind a square cylinder, we used the numerical basin that was 280 m long, 60 m wide and 5 m deep. Also a non-uniform grid containing 190 × 44 horizontal cells and 11 vertical layers was used as shown in Fig. 1. A constant flow propagates toward a square cylinder that is 30 m from the boundary inlet in which the drag coefficient of the square cylinder is 2.2 at a Reynolds number Re ~ 1.0 × 105 (Blevins, 1984). The constant velocity, V = (u, 0, 0) = 0.1 m/s, was specified at the boundary inlet. And the continuative boundary condition was applied at the flow outlet, while the smooth and the no-slip condition were used at the bottom boundary conditions. Here, all the calculation conditions were the same as in the previous study by Park et al. (2014).

FLOW-3D with the inclusion of the two different turbulence models, the k-ε and RNG k-ε models, was used to simulate wakes behind the square cylinder. Launder and Spalding (1974) and Yakhot et al. (1992) suggested the model constant Cμ = 0.09 (= Cμ0) and CμR = 0.085 (= CμR0) for the k-ε and RNG k-ε models, respectively. According to the experiment of Rodi (1972), the eddy viscosity for the simulation of wakes behind the cylinder must be 3 to 13 times higher than that for the round jet, plane jet and pipe flow. Thus we calculated wakes using constant values like 4Cμ0, 5 Cμ0, 6Cμ0, 7Cμ0, 8Cμ0 and 9Cμ0 or 4CμR0, 5CμR0, 6CμR0, 7 CμR0, 8CμR0 and 9CμR0 in both the turbulence two-equation models. In addition, the simulation was performed for over 5000 seconds with varying time steps controlled by the program. After the flow reached a steady state condition, the simulation results were analyzed (Park et al., 2014).

## 4. Model Results and Discussion

### 4.1 General Horizontal and Vertical Distribution of Wakes

Using the RNG k-ε model, we calculated the velocity fields, and the longitudinal (x-direction) velocity distribution is especially shown in Fig. 2. The wakes behind the square cylinder were similar to the long tail of a shooting star in the sky and the width of the 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 (Park et al., 2014).

From the longitudinal velocity along the centerline, we could find that a reverse flow occurred near the square cylinder due to the low pressure. After passing the cylinder, the flow rapidly recovered around x/ (CDD) = 10. The velocity gradient slowly varied till x/(CDD) = 50 and became stabilized after x/(CDD) = 50. These results were quite close to those obtained by Schlichting (1930). Resulting from the strong vertical mixing and interaction of the vortex, the differences in the horizontal velocity at the three layers were insignificant just behind the structure. These velocity characteristics looked very similar to the previous study by Park et al. (2014). On the other hand, the results of the k-ε model were similar to those of the RNG k-ε model, and so are not presented here.

### 4.2 The k-ε Model’s Results

The k-ε model, which solves two transport equations for the k and the ε, does not need the information of the mixing length. Although the eddy viscosity in the model depends on the k and ε, an empirical constant Cμ in eddy viscosity equation (Eq. 2) has to be specified. The standard Cμ (Cμ0 = 0.09) from the empirical observation was used everywhere (Launder and Spalding, 1974). Rodi (1972) noted that the Cμ for the wake experiments behind a sphere must be 3 to 13 times higher than that of for the jet or pipe flow experiments. Due to the fact that there were insufficient studies to suggest reasonable values of Cμ for the simulation of wakes behind a square cylinder, we conducted several scenarios using the k-ε model with different Cμ that ranged from 4Cμ0 to 9 Cμ0, found the best model constant and evaluated the applicability of the model constant for the simulation of wakes.

For all scenarios, the k-ε model results underestimated the wake width (Fig. 3) and the values of γ ranged from 0.843 to 0.877 (Table 1), which were worse than that of the one-equation model and mixing length model in the previous study by Park et al. (2014). Even if the EMVD was similar to that of the one-equation model, each velocity difference at three vertical layers was significant which meant that the model created weak vertical mixing. The RDC ranged from 56.2 % to 56.4 % in Table 1 which was relatively high when compared to the other models, so that the model had a limitation to require decent RDC at any Cμ . Despite the fact that the k-ε model employed two transport equations, the model could not successfully calculate wakes behind the square cylinder. The reason stemmed from the eddy viscosity hypothesis and ε-equation. The ε-equation included two empirical constants with suggested C1 = 1.44 and C2 = 1.92 in Eq. (5). It performed well for simple flows but not for complex flows as mentioned by Stephen (2000). To improve the accuracy of model calculation for wakes behind the square cylinder, the model constants must be adjusted and close investigations for their ranges are necessary.

### 4.3 The RNG k-ε Model’s Results

Many studies have contributed towards improving the performance of the k-ε model (Pope, 1978;Hanjalic and Launder, 1980;Bardina et al., 1983). Yakhot and Smith (1992) developed the RNG k-ε model based on theoretical grounds. The RNG k-ε model included an additional term in the ε-equation which could not be neglected and was important in the highly strained flows and low Reynolds number flows (Yakhot et al., 1992;Maurizi, 2000;Kim and Baik, 2004). Since the model is based on a theoretical ground, all the constants in the ε-equation derived from the theory itself and the values are C1 = 1.42 and C2 = 1.68. Therefore, the empirical constant CμR in the modified eddy viscosity equation plays a major role in calculating the eddy viscosity. The typical CμR (= CμR0) in the RNG k-ε model is 0.085, which is similar to the value Cμ0 (= 0.09) used in the k-ε model. We calculated the wakes behind the square cylinder using the RNG k-ε model with various ranges of CμR . The CμR ranges 4CμR0 to 9CμR0 were the same as those used in the k-ε model because of the closeness of two constant values as well as for the comparisons of two models.

As the CμR increased, the wake width increased and the increasing trend was similar to that of the results of the one-equation model in the previous study by Park et al. (2014). As shown Fig. 4, the range of γ was 0.963 to 0.995 (Table 2) which was the best when compared to other models. The EMVD showed large variations from -28.1 % to 22.4 % in Table 2 corresponding to the CμR . But the variations of EMVD become relatively insignificant like 5 % to 7 % when the CμR changed from 6CμR0 to 8CμR0. At these ranges, the model predictions agreed well with the analytic solutions. The RDC generally increased along with the increase of CμR and its range was 0.662 to 1.315 (Table 2). The best fit of the RDC was 1.032 which showed 3.2 % error when it occurred at CμR = 6CμR0 as shown in Fig. 4 and Table 2. The AME (Absolute Mean Error) with the γ, EMVD and RDC ranged from 3.8 % to 19.3 % (Table 2) and these values were much better than those of the k-ε model (Fig. 3 and Fig. 4). This confirmed that the additional term and the model constants derived from the theory in the ε-equation were key components in simulating the wakes behind the square cylinder. From the RNG k-ε model results, the agreement of model predictions and analytic solutions was exceptionally good when the CμR was 6 times higher than the suggested value. In addition, on the basis of the experiment of Rodi (1972) at P/ε < 0.5, the eddy viscosity for wakes behind the cylinder must be 3 to 13 times higher than that for the round jet, plane jet and pipe flow simulations. As a result, the high value of CμR might be appropriate for the simulation of wakes behind the square cylinder.

## 5. Conclusions

The study simulated wakes behind a square cylinder using two-equation turbulence models based on the eddy viscosity concept, the k-ε and RNG k-ε models, and investigated the applicable range of the major empirical constant in each model. For the comparisons between the model predictions and analytical solutions, we employed three skill assessments:, the correlation coefficient for the similarity of the wake shape, the error of maximum velocity difference for the wake velocity distribution, and the ratio of drag coefficient for the flow patterns. Based on the above results, the following conclusions were arrived at.

The k-ε model results underestimated the wake width, and its performance was worse than the one-equation and mixing length models in the author’s previous study. The reason stemmed from the eddy viscosity hypothesis and ε-equation. The use of the ε -equation, including empirical constants with suggested values, might be inappropriate for the simulation of wakes behind the square cylinder. Therefore, adjusting the model constants was necessary for the accuracy of the model calculations.

Excellent results were obtained with the RNG k-ε model, resulting from an additional term and the model constants derived from the theory in the ε-equation. The AME averaging the γ, EMVD and RDC were less than 3.8 % at CμR = 6CμR0, which was 6 times higher than the suggested value. It indicated that the high value of CμR might be appropriate for the simulation of wakes behind the square cylinder.

## Acknowledgements

This study was financially supported by Chonnam National University, 2011.

## Figure

Diagram of a square cylinder, a water basin and grid systems (after Park et al., 2014).

Velocity distributions on the surface (k = 11), middle (k = 6) and bottom (k = 1) layer for the RNG k-ε model where k is the vertical layer grid number.

Comparison between the analytical and numerical solutions by the k-ε model.

Comparison between the analytical and numerical solutions by the RNG k-ε model.

## Table

γ, EMDV, RDC and AME from the k-ε model’s results (Cμ0 = 0.09)

γ, EMDV, RDC and AME from the RNG k-ε model’s results (CμR0 = 0.085)

## Reference

1. Bardina, J. , Ferziger, J. H. and Reynolds, W. C(1983), Improved Turbulence Models Based on Large Eddy Simulation of Homogeneous, Incompressible, Turbulent Flows, Thermosciences Division, Department of Mechanical Engineering, Report No. TF-19, Stanford University, USA, p. 193.
2. Blevins, R. D. (1984), Applied Fluid Dynamics Handbook, van Nostrand Reinhold Company, pp. 279-311.
3. Bosch, G. and Rodi, W. (1998), Simulation of Vortex Shedding Past a Square Cylinder with Different Turbulence Models, International Journal for Numerical Methods in Fluids, Vol. 28, pp. 601-616.
4. Chen, L. , J. Y. Tub and G. H. Yeohc(2003), Numerical Simulation of Turbulent Wake Flows behind Two Side-by-Side Cylinders, Journal of Fluids and Structures, Vol. 18, pp. 387-403.
5. Chou, P. Y. (1945), On Velocity Correlations and the Solutions of the Equation of Turbulent Fluctuation, Applied Mathematics, Vol. 111, No. 1, pp. 38-54.
6. Ettema, R. , B. Yoon, T. Nakato and M. Muste(2004), A Review of Scour Conditions and Scour-Estimation Difficulties for Bridge Abutments, KSCE Journal of Civil Engineering, Vol. 8, No. 6, pp. 643-650.
7. Flow Science(1999), FLOW-3D User's Manual, Chapter 3 Theory, Los Alamos, NM, USA, pp. 1-164.
8. Görtler, H. (1942), Berechnung von Aufgaben der Freien Turbulenz auf Grund Eines Neuen Näherungsa-Ansatzes, ZAMM, Vol. 22, pp. 244-245.
9. Hanjalic, K. and B. E. Launder(1980), Sensitizing the Dissipation Equation to Irrotational Strains, J. Fluids Eng. Trans. ASME, Vol. 102, pp. 34-40.
10. Honzejk, V. and K. Frana(2008), A Turbulent Flow Past a cylinder, Journal of Applied Science in the Thermo Dynamics and Fluid Mechanics, Vol. 2, No. 22, pp. 1-6.
11. Ishino, K. , H. Otani, R. Okada and Y. Nakagawa(1993), The Flow Structure around a Cylindrical Pier for the Flow of Transcritical Reynolds Number, Proc. of XXV Congress of IAHR, Vol. 5, pp. 417-424.
12. Kim, J. J. and J. J. Baik(2004), A Numerical Study of the Effects of Ambient Wind Direction on Flow and Dispersion in Urban Street Canyons Using the RNG k-ε Turbulence Model. Atmosphere Environment, Vol. 38, pp. 3039-3048.
13. Launder, B. E. and D. B. Spalding(1974), The Numerical Computation of Turbulent Flows, Computer Methods in Applied mechanics and Engineering, Vol. 3, pp. 269-289.
14. Liou, T. M. , S. H. Chen and P. W. Hwang(2002), Large Eddy Simulation of Turbulent Wake behind a Square Cylinder with a Nearby Wall, Journal of Fluids Engineering, Vol. 124, pp. 81-90.
15. Maurizi, A. (2000), Numerical Simulation of Turbulent Flows over 2-D Valleys Using Three Versions of the k-ε Closure Model, Journal of Wind Engineering and Industrial Aerodynamics, Vol. 85, pp. 59-73.
16. Olsen, N. R. B. and M. C. Melaaen(1993), Three-Dimensional Calculation of Scour around Cylinders, Journal of Hydraulic Engineering, ASCE, Vol. 119, No. 9, pp. 1048-1054.
17. Park, I. H. , Y. J. Cho, T. Y. Kim, M. O. Lee and S. S. Hwang(2014), Analysis of Empirical Constant of Eddy Viscosity by Zero- and One-Equation Model in Wake Simulation. Journal of the Korean Society of Marine Environment & Safety, Vol. 20, No. 3, pp. 323-333.
18. Park, I. H. , G. H. Lee and Y. J. Cho(2009), Drag Coefficient Estimation of Pile Type Structures by Numerical Water Basin Experiments, Journal of Korean Society of Coastal and Ocean Engineering, Vol. 21, No. 1, pp. 45-53.
19. Park, I. H. , J. S. Lee and M. O. Lee(1998), A Numerical Model of Large Scale Grid for Two-Dimensional Wake behind Bodies, Journal of Korean Society of Coastal and Ocean Engineering, Vol. 10, No. 2, pp. 83-92.
20. Pope, S. B. (1978), An Explanation of the Turbulent Round-Jet/Plane-Jet Anomaly, AIAA J., Vol. 16, pp. 279-281.
21. Reichardt, H. (1951), Gesetzmässigkeiten der Freien Turbulenz, VDI-Forschung-Sheft, 414 (1942), 2nd Ed., p. 1951.
22. Richardson, J. E. and V. G. Panchang(1998), Three- Dimensional Simulation of Scour-Inducing Flow at Bridge Pipers, Journal of Hydraulic Engineering, ASCE, Vol. 124, No. 9, pp. 530-540.
23. Rodi, W. (1972), The Prediction of Free Turbulent Boundary Layers by Use of a Two-Equation Model of Turbulence, Thesis Submitted for the Degree of Doctor of Philosophy in the Faculty of Engineering, University of London and for the Diploma of Membership of Imperial College, p. 310.
24. Schlichting, H. (1930), Über das Ebene Windschatten Problem. Diss Göttingen, Ing.-Arch., Vol. 1, pp. 533-571.
25. Schlichting, H. (1979), Boundary Layer Theory, 7th Ed., McGraw-Hill, New York (Translated by J. Kestin), p. 817.
26. Stephen, B. P. (2000), Turbulent Flows, Cambridge University Press, p. 771.
27. Tseng, M. H. , C. L. Yen and C. C. S. Song(2000), Computation of Three-Dimensional Flow around Square and Circular Piers, International Journal for Numerical Methods in Fluids, Vol. 34, pp. 207-227.
28. Yakhot, V. and S. A. Orszag(1986), Renormalization-group analysis of turbulence, Journal of Scientific Computing, Vol. 1, pp. 3-51.
29. Yakhot, V. , S. A. Orszag, S. Thangam, T. B. Gatski and C. G. Speziale(1992), Development of Turbulence Models for Shear Flows by a Double Expansion Technique, Phys. Fluid. A, Vol. 4, No. 7, pp. 1510-1520.
30. Yakhot, V. and L. M. Smith(1992), The Re-Normalization Group, the ε-Expansion and Derivation of Turbulence Models, Journal of Scientific Computing, Vol. 7, pp. 35-61.
31. Yang, W. and S. U. Choi(2002), Three-Dimensional Numerical Simulation of Local Scour around the Bridge Pier Using Large Eddy Simulation. Journal of Korean Society of Civil Engineering, Vol. 22 (4-B), pp. 437-446.