Numerical simulation of flow past a square cylinder using Partially-Averaged Navier–Stokes model

ARTICLE IN PRESS Journal of Wind Engineering and Industrial Aerodynamics 97 (2009) 37–47

Contents lists available at ScienceDirect

Journal of Wind Engineering and Industrial Aerodynamics journal homepage: www.elsevier.com/locate/jweia

Numerical simulation of flow past a square cylinder using Partially-Averaged Navier–Stokes model Chi-Su Song, Seung-O Park  Department of Aerospace Engineering, Korea Advanced Institute of Science and Technology, 373-1, Guseong-dong Yuseong-gu, Daejeon 305-701, Korea

a r t i c l e in fo

abstract

Article history: Received 8 November 2007 Received in revised form 12 November 2008 Accepted 14 November 2008 Available online 6 January 2009

In this study, we carry out a two-stage partially averaged Navier–Stokes (PANS) simulation past a square cylinder to examine the effect of a grid system on the predicted results. The important focus in PANS is how the control parameter which is the ratio of unresolved-to-total turbulent kinetic energy is determined in the flow field. In this work, the parameter is evaluated by using an equation that we propose. The model equation utilizes the turbulent length scale, the Kolmogorov length scale and the cut-off grid length. We first perform RANS simulation to approximately evaluate the turbulent length scale and the Kolmogorov length scale. Then we construct five different grid systems based on these length scales and perform PANS simulation in which the length scales to determine control parameter are also updated as the simulation proceeds. Simulation results are compared with experimental data and the other LES and DES results. The comparison suggests that PANS approach is very effective in simulating separated turbulent flows in the sense that accurate predictions can be obtained by using a much coarser grid system than is required for LES. We then suggest the grid spacing requirements for a proper PANS simulation. & 2008 Elsevier Ltd. All rights reserved.

Keywords: PANS Square cylinder Control parameter Two-stage approach Cut-off grid length

1. Introduction Flow past a square cylinder has been studied by many investigators for a long period of time. As is well known, flow separation and vortex shedding are typical features in this particular flow. Until recently, numerical investigations of this square cylinder flow have adopted RANS (Reynolds Averaged Navier–Stokes) equations with turbulence model. The k–e turbulence model has been most widely used (Bosch and Rodi, 1998; Lee, 1997). However, this approach has often predicted results that are inconsistent with experimental data (Lyn et al., 1995). LES (Large Eddy Simulation) is of course a much better approach for calculation of unsteady turbulent flows. In LES, only small scale motion is modeled while large scale eddies are directly resolved. The flow past a square cylinder is often exemplified as a benchmark problem for testing various sub-grid scale models of LES (Rodi et al., 1997; Sohankar et al., 2000; Voke, 1997). However, because of the overwhelming requirement of computing cost owing to many grid points required especially near the wall, LES is not yet a practical tool for engineering application. Recently, an alternative method called Hybrid RANS/LES has been proposed to efficiently combine the advantages of RANS and LES. Spalart et al.

 Corresponding author. Tel.: +82 42 869 3713; fax: +82 42 869 3710.

E-mail address: [email protected] (S.-O. Park). 0167-6105/$ - see front matter & 2008 Elsevier Ltd. All rights reserved. doi:10.1016/j.jweia.2008.11.004

(1997) proposed DES (Detached Eddy Simulation) in which the length scale of Spalart–Allmaras turbulence model (Spalart and Allmaras, 1992) is determined as a minimum value between wall distance and grid size so that the model can calculate eddy viscosity near the wall and the sub-grid scale eddy viscosity away from the wall. It is well accepted that DES performs well for flows with massive separation. However, Nikitin et al. (2000) in their numerical study of channel flow with DES pointed out that mean velocity profiles had mismatch between RANS and LES region; the velocity gradients were too steep near interface of RANS and LES region. Speziale (1998) proposed a VLES (Very Large Eddy Simulation) method which allows computation smoothly bridging from RANS simulation to DNS (Direct Numerical Simulation) depending on the grid resolution. But this approach requires a substantially more sophisticated modeling approach such as adjusting model constants (Fro¨hlich and Terzi, 2008). Girimaji et al. (2003) suggested Partially Averaged Navier–Stokes (PANS) method based on the ratio of the unresolved to total kinetic energy (k) and dissipation rate (e). Recently, Fro¨hlich and Terzi (2008) emphasized that a main advantage of the PANS model is its easy implementation into an existing RANS solver as only the turbulence model coefficients need to be changed. The PANS model was tested against some typical turbulent flows such as the flow around square cylinder and nozzle flows (Abdol-Hamid and Elmiligui, 2005; Girimaji, 2006; Girimaji et al., 2006). The predictions were found to be challengeable. Lakshmipathy and

ARTICLE IN PRESS 38

C.-S. Song, S.-O. Park / J. Wind Eng. Ind. Aerodyn. 97 (2009) 37–47

Girimaji (2006) further derived PANS equations based on the k–o model, which was used to simulate the flow past a circular cylinder and a backward facing step flow in their work. In this study, we carry out a PANS simulation past a square cylinder to examine the effect of grid system on the predicted results. For a hybrid RANS/LES simulation in general, grid dependency of the solution poses a serious problem for practitioners. As the simulation is dictated by grid system, we often pose equations such as how good the solution is and what would be the minimum requirement for grid spacing and so on. These questions become even more serious when we need to deal with large scale flows such as building aerodynamics. Issues concerning grid spacing have been raised by many investigators so far. However, a systematic investigation seems to be rare. In this work, we construct various grid systems based on the length scale obtained from a RANS solution with two equation turbulence model. The results of the simulation are then compared with experimental data, another LES (Lubcke et al., 2001; Sohankar et al., 2000) or DES (Barone and Roy, 2006; Schmidt and Thiele, 2002) results. The advantage of PANS procedure compared to LES is then elucidated in the view point of practical engineering calculation.

2. PANS model 2.1. PANS equations In this section, we briefly summarize PANS equations for the convenience of presentation. Girimaji et al. (2003) first derived PANS equations based on k–e model. Later, Lakshmipathy and Girimaji (2006) suggested PANS equations based on k–o model. They showed that the k–o based PANS performed better than the k–e model based one for the flow past a two-dimensional circular cylinder. As the flow past a square cylinder is similar in character to the flow past a two-dimensional circular cylinder, we choose to adopt PANS model based on k–o model. The set of governing equations are:

qU i ¼0 qxi

(1)

qU i qU i 1 qP q þ Uj ¼ þ ½2nS ij  tij  qt qxj r qxi qxj

(2)

qku qku q  þ Uj ¼ P u  b k u ou þ qt qxj qxj



nþ



nt qku sku qxj



(3)

The model constant b^ in Eq. (4) is given by   1  b^ ¼ ab þ ðb  ab Þ fo where fo ¼ fe/fk. The control parameters fk and fe are the following ratios: f k ¼ ku =k;

f  ¼ u =

(9)

The k and e are the total turbulent kinetic energy and dissipation rate, respectively. In earlier works (Girimaji et al., 2003, 2006; Girimaji, 2006; Lakshmipathy and Girimaji, 2006), fk was often dealt with as a constant between 0 and 1, whereas fe was taken to be unity. With fe ¼ 1, we simply have 1/fo ¼ fk and ou ¼ o/fk.

2.2. Control parameter fk The control parameter, fk, which is the ratio of the unresolved (or modeled) turbulent kinetic energy to the total turbulent kinetic energy, is decisive in PANS simulation as can be seen in the set of governing equations. Unless we have full details of turbulence for a specific flow, we can hardly determine fk precisely, which is not possible in practice. Thus, we have either to determine fk prior to the simulation or to update it as the simulation proceeds. In the early stage of PANS application, fk was assumed to be constant such as 0.7 or 0.4, etc (Girimaji et al., 2003, 2006; Girimaji, 2006; Lakshmipathy and Girimaji, 2006). As the assumption of constant fk is not reasonable, some formula has been suggested to take into account of the grid size and turbulence length scale. Girimaji and Abdol-Hamid (2005) proposed the following form to determine fk   D 2=3 fk  3 lturb

(10)

where lturb ¼ k1.5/e, and D ¼ min(Dx, Dy, Dz). Elmiligui et al. (2004) introduced a variable form of fk in space. Abdol-Hamid and Girimaji (2004) suggested a ‘‘two-stage approach’’ to distribute the control parameter for PANS simulation. In this work, we adopt a two-stage process to determine fk according to our proposal as described below. We assume that the grid spacing for PANS simulation is located between the turbulent length scale (lturb) and the Kolmogorov length scale (Z) as shown in Fig. 1. By assuming that all the length scales appearing in Fig. 1 belong to the inertial range, we have,

0 1 Rk  2=3    2=3  2=3  Z 2=3  Z 2=3  kZD EðkÞ k2=3  kD D 2=3 B 1  DZ D Z C ¼ 2=3 f k ¼ R kl ¼ 1  ’ @ A  2=3 lturb Z D lturb  kZ EðkÞ kZ  kl2=3 1 Z

(11)

lturb



qou qou ou q qou þ Uj ¼a P u  b^ o2u þ ðn þ su nt Þ qt qxj ku qxj qxj



(4)

The subscript u refers to the unresolved quantities. The Pu in Eqs. (3) and (4) is the unresolved scale production term given by P u ¼ tij S ij

(5)

The unresolved Reynolds stress (tij) is modeled as:

tij ¼ 2nt S ij  23ku dij nt ¼

(6)

ku

(7)

ou

The model constants in the equations above are:

b ¼ 0:09;

a ¼ 5=9;

b ¼ 0:075;

sku ¼ 0:5;

su ¼ 2:0.

(8)

Fig. 1. Energy spectrum.

ARTICLE IN PRESS C.-S. Song, S.-O. Park / J. Wind Eng. Ind. Aerodyn. 97 (2009) 37–47

where D is the cut-off grid length often defined as D ¼ max(Dx, Dy, Dz) (Spalart et al., 1997; Schmidt and Thiele, 2002) or D ¼ (DxDyDz)1/3 (Rodi, 1997; Sohankar et al., 2000). To determine fk, we need to estimate lturb and Z as functions of space coordinates. For this purpose, we first perform RANS calculation to determine k and o. With estimates of lturb and Z available, we construct a suitable grid system to carry out PANS simulation. In the next section, we present how we apply this procedure to simulate the flow past a square cylinder. When the control parameter (fk) varies in time and space, we need to take this into consideration in the transport equation of ku and ou. It can be readily shown that Eqs. (3) and (4) should be changed to      qku qku q nu qku qf k qf  þk ¼ P u  b ku ou þ þn þ Uj þ Uj k qxj sku qt qxj qxj qt qxj

39

largest and the Kolmogorov length scale is smallest in the wake region. The spatially averaged value of the normalized Kolmogorov length scale (Zavg/D) over 5pz/Dp5 and 1px/D is about 0.003 for this particular flow. For PANS simulation, we select four reference values of D/D to design grid systems. They are: 0.3 (D ¼ 100Zavg), 0.15 (D ¼ 50Zavg), 0.1 (D ¼ 33Zavg), and 0.075 (D ¼ 25Zavg). The cutoff grid length (D) is chosen according to D ¼ max(Dx, Dy, Dz). Therefore D works as an upper limit for grid spacing in each direction. In devising a grid system for turbulent flow simulations, a focus region in computational domain where turbulent action is most significant need be identified for the economy of computation.

(12)     qou qou ou q qou qf o qf þ Uj þ Uj o þo ¼a P u  b^ o2u þ ðn þ su nt Þ qt qxj ku qxj qt qxj qxj

(13) In terms of fk, ku and ou, the length scales in Eq. (11) are expressed as 

lturb ¼ k1=2 u =ðb ou Þ

(14)

Z ¼ ðn3 =Þ1=4 ¼ ðn3 =ðb ou ku ÞÞ1=4

(15)

Fig. 2. Computational domain for RANS simulation.

In the present ‘‘two-stage’’ process, initial distribution of fk is determined from the steady RANS solution. And then Eqs. (1), (2), (12), and (13) are integrated and fk distribution is updated by using Eqs. (11), (14), and (15).

3. Flow past a square cylinder For the PANS simulation of the flow past a square cylinder, we first perform steady RANS simulation with the k–o turbulence model. Various grid arrangements for this calculation as given in Table 1 are tested. The computational domain is shown in Fig. 2. Fig. 3 displays the distributions of the mean velocity U along the center line together with the experimental data from Lyn et al. (1995). Based on the results of Fig. 3, the grid system ‘‘st-3’’ is judged to be good enough for RANS simulation. To devise a suitable distribution of grid spacing for PANS simulation, we need to estimate distributions of lturb and Z from the RANS simulation results, which are illustrated in Fig. 4. Fig. 4 indicates that the normalized Kolmogorov length scale (Z/D) ranges from 103 to 102 while the normalized turbulent length scale (lturb/D) is in the range of 101–100. We also see that the turbulent length scale is

Fig. 3. U velocity distribution along the center line.

Table 1 Grid system for two-dimensional RANS simulation. Cases

Total grid (Nx  Nz)

Cylinder’s resolution (Nx  Nz)

Maximum and minimum resolution

Stretching ratio qx ¼ Dxi/Dxi1 qz ¼ Dzi/Dzi1

Dxmin ¼ Dzmin ¼ 0.0065D, Dxmax ¼ 2.5D, Dzmax ¼ 1.3D Dxmin ¼ Dzmin ¼ 0.0052D, Dxmax ¼ 1.9D, Dzmax ¼ 0.9D Dxmin ¼ Dzmin ¼ 0.0048D, Dxmax ¼ 1.4D, Dzmax ¼ 0.75D Dxmin ¼ Dzmin ¼ 0.004D, Dxmax ¼ 1.2D, Dzmax ¼ 0.6D Dxmin ¼ Dzmin ¼ 0.0035D, Dxmax ¼ 0.9D, Dzmax ¼ 0.5D

Max Max Max Max Max Max Max Max Max Max

Case st-1

65  52

17  17

Case st-2

91  72

23  23

Case st-3

109  86

27  27

Case st-4

127  102

31  31

Case st-5

154  122

35  35

qx ¼ 1.2 qz ¼ 1.48 qx ¼ 1.14 qz ¼ 1.3 qx ¼ 1.13 qz ¼ 1.25 qx ¼ 1.1 qz ¼ 1.21 qx ¼ 1.09 qz ¼ 1.15

ARTICLE IN PRESS 40

C.-S. Song, S.-O. Park / J. Wind Eng. Ind. Aerodyn. 97 (2009) 37–47

Fig. 4. The distribution of turbulent length scale (lturb/D; solid symbol) and Kolmogorov length scale (Z/D; open symbol) at several stations, D being the width of the square cylinder.

normally taken to be uniform, is usually larger than the grid spacing in the other directions. For the spanwise length of 4D, the four preselected D’s give the number of grid points to be 14, 27, 40, and 54, as can be seen in Table 2. The turbulent length scale, Kolmogorov length scale, and grid spacing for the present PANS simulation cases in streamwise direction (x) and vertical direction (z) are illustrated in Fig. 6. We see that all the grid spacings corresponding to the cases of Table 2 are located between the turbulent length scale and the Kolmogorov length scale in the computational region of our interest. To examine the effect of spanwise grid spacing on the simulation, case A1 was added as given in Table 2. The spanswise grid spacing of case A1 is about half of that of case A2. The cut-off grid length (D/D) in this case is 0.04 for x/Do3.0. From Fig. 4, we find that the turbulent length scale is approximately in the range of 102Z–103Z. By using Eq. (11), we can easily construct the range of fk as a function of D/Zavg as shown in Fig. 7, which also illustrates the range of fk values for various cases of grid system selected in the present work. For case A1, the range shown in the figure is for the region x/Do3 corresponding to D/Zavg ¼ 13.3. The initial flow field for the simulation is obtained from the preliminary RANS simulation results. Initial fields of ku and ou are calculated from k and o values of the RANS result by using fk distributions estimated by Eq. (11) for each grid system. Periodic boundary condition is employed in the spanwise direction. Noslip boundary condition for velocity and ku ¼ 0, wu ¼ 60n/(b(y1)2) for turbulence quantities are imposed at solid boundary; y1 is the physical distance of the first point away from the wall (Baldina et al., 1997). The Neumann boundary condition is adopted at the upper and lower boundaries. At the inlet boundary, uniform flow with 2% turbulence intensity in accordance with the experimental condition (Lyn et al., 1995) is specified. We choose the nondimensional time step, defined by tU/D, of 0.005, which is found to be sufficient to track the unsteady characteristics of the flow. The three-dimensional incompressible Navier–Stokes equations are integrated using a fully implicit fractional step method. The governing equations are transformed to generalized curvilinear coordinates on a non-staggered grid. The convection term is discretized using a fifth order upwind scheme and all the other terms in momentum and pressure-Poisson equation are approximated using the second-order central difference scheme. Time integration is done using dual-time-stepping algorithm and local time stepping method is adopted (Constantinescu and Squires, 2004).

4. Results and discussion Fig. 5. Computational domain with focus region denoted by dashed lines.

In the present work, we select this focus region as the region confined by 0ox/Do7 and 2oz/Do2 after examining the experimental data of Lyn et al. (1995). The computational domain including the focus region is sketched in Fig. 5. As PANS simulation requires three-dimensional domain, we select the spanwise dimension to be 4D which is a typical spanwise length widely chosen in the other LES or DES simulations (Barone and Roy, 2006; Lubcke et al., 2001; Rodi, 1997; Rodi et al., 1997; Schmidt and Thiele, 2002; Sohankar et al., 2000; Voke, 1997). Details of various grid systems adopted in the present work are summarized in Table 2. The Dzþ s expressed in wall unit in the min table are based on the wall shear stress at the center point of the upper surface of the square cylinder (x/D ¼ z/D ¼ 0.5). When the mean flow is two-dimensional, the spanwise grid spacing,

The simulation for all the cases is carried out for 260 nondimensional time period which covers about 35 vortex shedding cycles. It is found that this period is sufficient to get a converged time-averaged flow field. Table 3 summarizes typical flow parameters from the present simulation including the recirculation length (lR/D), the Strouhal number (St) of the vortex shedding, the mean drag coefficient (Cd), and the rms of Cd fluctuation. Comparisons are made with the experimental data of Lyn et al. (1995), LES results of Sohankar et al. (2000). Also included in the table are LES data such as UMIST-2, UKAHY-2, and TAMU-2 as compiled by Rodi et al. (1997) and NT7, UOI, TIT, GRO as collected by Voke (1997), DES result (DES-fine) of Barone and Roy (2006) and Schmidt and Thiele (2002). An important aspect to be kept in mind when comparing predicted results of various simulations is about the grid system employed. For all of the cases listed in Table 3 except for the case of ‘LES-OEDSMF’, the cross-sections of the computational domain are the same: the spanwise length is 4D

ARTICLE IN PRESS C.-S. Song, S.-O. Park / J. Wind Eng. Ind. Aerodyn. 97 (2009) 37–47

41

Table 2 Grid system for PANS simulation. Cases

Grid size of computational domain and of focus region (Nx  Nz)

Cylinder’s resolution (Nx  Nz)

Spanwise grid resolution

Minimum and maximum resolution in focus region

Upper bound of cut-off grid size in the focus region

Stretching ratio qx ¼ Dxi/Dxi1 qz ¼ Dzi/Dzi1

Case A1

204  122 and 160  100

40  40

Ny ¼ 100

Dxmin ¼ Dzmin ¼ 0.003D, Dzþ min ¼ 0:9

D ¼ 25Zavg ( ¼ 0.075D)

Max qx ¼ 1.03

Dy ¼ 0.04D Ny ¼ 54

Case A2

Max qz ¼ 1.13

Dxmax ¼ 0.073D, Dzmax ¼ 0.074D

DyE0.07D Case B

187  114 and 134  95

35  35

Ny ¼ 40

Dy ¼ 0.1D Case C

174  108 and 124  89

30  30

Ny ¼ 27

DyE0.15D Case D

146  86 and 86  67

27  27

Ny ¼ 14

DyE0.29D RANS(2D) 109  86

27  27

–

Dxmin ¼ Dzmin ¼ 0.0035D, Dzþ ¼ 1:0 min Dxmax ¼ 0.1D, Dzmax ¼ 0.08D

D ¼ 33Zavg ( ¼ 0.10D)

Dxmin ¼ Dzmin ¼ 0.004D, Dzþ ¼ 1:2 min Dxmax ¼ 0.15D, Dzmax ¼ 0.1D

D ¼ 50Zavg ( ¼ 0.15D)

Max qz ¼ 1.14

Max qx ¼ 1.06 Max qz ¼ 1.19

D ¼ 100Zavg ( ¼ 0.3D) Dxmin ¼ Dzmin ¼ 0.0048D, Dzþ ¼ 1:5 min Dxmax ¼ 0.28D, Dzmax ¼ 0.29D Dxmin ¼ Dzmin ¼ 0.0048D, Dxmax ¼ 1.4D, Dzmax ¼ 0.75D

Max qx ¼ 1.04

–

Max qx ¼ 1.09 Max qz ¼ 1.27

–

Fig. 6. Grid spacing, turbulent length scale, and Kolmogorov length scale distributions in streamwise and vertical directions in the focus region: (a) streamwise direction and (b) vertical direction.

and the height is 14D (see Fig. 5). For the case of LES-OEDSMF, the cross-section is 4D  15.7D. The streamwise length of the computational domain varies case by case as can be seen from Table 3. The lengths Lu and Ld in the table refer to the upstream length between the inlet boundary and the front surface of the cylinder and the downstream length between the rear surface of cylinder and the outlet boundary, respectively. We see that the computational domain of the present cases is largest among the other cases listed except for the cases of ‘DES-fine’. We also find that the grid arrangements for the present simulations are much coarser than that of ‘DES-fine’. The grid systems of ‘Case B’ and ‘Case C’ are comparable to that of ‘DES-A’. It is worthwhile to note that the LES simulation cases with wall function adopted much

coarser grid system than the other cases contained in the table. It should be pointed out that Kogaki et al. (1997) recommended to provide 40 node points (10 node points per cylinder length) and Voke (1997) suggested to have at least 32 meshes to adequately cover the spanwise length of 4 D for this particular flow. In this regard, case A1 and A2, and B of the present work satisfy this suggestion. From Table 3, we find that the length of the recirculation bubble (lR/D) varies most widely with different choices of the grid system while the drag coefficient, Cd is least sensitive to the grid system. The mean drag coefficients (Cd) and the Strouhal number (St) from all of the present PANS simulations agree well with the experimental data. This feature is in accordance with Rodi’s

ARTICLE IN PRESS 42

C.-S. Song, S.-O. Park / J. Wind Eng. Ind. Aerodyn. 97 (2009) 37–47

observation (Rodi et al., 1997; Rodi, 1997) that these values are relatively insensitive to the grid system. The distributions of the time and spanwise averaged streamwise velocity U along the centerline are displayed in Fig. 8. LES result (LES-OEDSMF) (Sohankar et al., 2000) included in the figure show that their simulation overpredicts much the velocity which is associated with a rather short recirculation length (Table 3). Whereas, DES result (DES-fine) of Barone and Roy (2006) is seen to be in much better agreement with the experimental data. The length of recirculation and the streamwise velocity profiles including the maximum negative velocity are far better predicted in this DES than LES (Sohankar et al., 2000) as given in Table 3. For the present cases of PANS simulation, the data from the cases A1, A2, and B are judged to be in good agreement with the experimental data. The cut-off grid length (D) of case A1 is much smaller than that of case A2 in the region x/Do3. A careful look at

Fig. 7. The range of fk (dashed line) with lturb variations for various grid systems.

Fig. 8 and the data of Table 3 suggests that the maximum negative velocity in the recirculation region gets larger as the spanwise grid spacing becomes smaller. The maximum negative velocity of case A1 agrees best with the experimental data. In spite of the excellent agreement with the experimental data illustrated in Fig. 8, we find however that the recirculation length of case A1 agrees with the experimental data no better than that of case A2. From the data of Table 3, we may easily notice that the recirculation length is not at all correlated consistently with the spanwise grid spacing. We also find from Fig. 8 that the case D is not good enough in the sense that the negative peak velocity in the recirculation bubble as well as the recirculation length (Table 3) and the recovery tendency of the streamwise velocity profile in the wake region are not in favorable agreement with the experimental data when compared to the other simulation cases of the present work. For these reasons stated above, we exclude case D from further discussion. We note, however, that the results of case D are as good as or better than those of the cases of ‘LES with wall function’. We compare the time and spanwise averaged longitudinal velocity profiles at several downstream locations (x/D ¼ 0.5, 2.5, 4.5, and 6.5) with the experimental data of Lyn et al. (1995) and LES data of Lubcke et al. (2001) in Fig. 9. We find that the predicted velocity profiles shown are in good agreement with the experimental data. The rms velocity fluctuations along the center line are presented in Fig. 10. It is seen that the present simulation data are in good agreement with the experimental data. The resolved rms velocity fluctuations increase slightly as the grid becomes finer as expected. The averaged pressure coefficient (CP) distribution around the square cylinder with the experimental data is presented in Fig. 11. We choose the pressure at the point of x/D ¼ 4.5 and z/D ¼ 0 in the free stream region as the reference free stream pressure to estimate CP. Fig. 11 indicates that the pressure coefficients on the front side facing incoming flow are fairly well predicted in all the cases of the present PANS simulation, whereas the pressure distributions on the rear side are slightly overpredicted than the experimental data of Bearman and Obasaju (1982). The pressure distribution on the rear side (face CD) and maximum negative

Table 3 Summary of typical parameters from various studies. Authors

Method

Case

St.

Cd

lR/D

rms Cd

(Nx, Ny, Nz) or Nx  Nz and Ny

Computational domain size (Lu/D, Ld/D)

B.C. at wall

Rodi et al. (1997)

LES

UKAHY2

0.130

2.30

0.96

0.14

(146, 20, 146)

(4.5, 13.5)

UMIST2

0.090

2.02

0.71

–

(140, 13, 81)

(4.5, 14.5)

TAMU2

0.140

2.77

0.44

0.19

(165, 17, 113)

(4.5, 14.5)

NT7

0.131

2.05

0.89

0.12

(140, 32, 103)

(4.5, 14.0)

UOI TIT GRO LESOEDSMF DES-fine

0.130 0.131 0.133 0.132

2.03 2.62 2.09 2.09

0.70 0.73 1.11 0.58

0.18 0.23 0.18 0.20

(192, 48, 160) (121, 127, 113) (280, 64, 210) (265, 25, 161)

(4.5, 13.5) (4.5, 9.5) (4.5, 19.5) (7.4, 15.8)

Wall function Wall function Wall function Wall function No-slip No-slip No-slip No-slip

0.125

2.11

0.82

0.26

(10.0, 20.0)

No-slip

DES-A

0.130

2.42

0.66

0.28

(5.0, 15.0)

No-slip

Case Case Case Case Case –

0.133 0.133 0.133 0.134 0.136 0.132

2.09 2.08 2.08 2.06 1.99 2.10

0.83 0.85 0.94 0.96 0.79 0.88

0.23 0.23 0.22 0.20 0.18 –

Nx  Nz ¼ 88, 200, Ny ¼ 96 Nx  Nz ¼ 32, 300, Ny ¼ 20 (204, 100, 122) (204, 54, 122) (187, 40, 114) (174, 27, 108) (146, 14, 86) –

(5.0, (5.0, (5.0, (5.0, (5.0, –

No-slip No-slip No-slip No-slip No-slip –

Voke (1997)

Sohankar et al. (2000) Barone and Roy (2006) Schmidt and Thiele (2002) Present

Lyn et al. (1995)

DES

PANS

Exp.

A1 A2 B C D

20.0) 20.0) 20.0) 20.0) 20.0)

ARTICLE IN PRESS C.-S. Song, S.-O. Park / J. Wind Eng. Ind. Aerodyn. 97 (2009) 37–47

43

Fig. 8. Streamwise velocity (U) distribution along the center line: (a) wake region and (b) recirculation region.

Fig. 9. Streamwise velocity profiles: (a) x/D ¼ 0.5, (b) x/D ¼ 2.5, (c) x/D ¼ 4.5, and (d) x/D ¼ 6.5.

velocity (see Fig. 8) predicted by using the grid system A1 is in excellent agreement with the experimental data, which eventually results in a better prediction of the mean drag coefficient

(Table 3). We comment here again that the finer spanwise grid spacing of case A1 makes the difference when compared to case A2. It is worthwhile to comment here about the length of time

ARTICLE IN PRESS 44

C.-S. Song, S.-O. Park / J. Wind Eng. Ind. Aerodyn. 97 (2009) 37–47

Fig. 10. Predicted normal stress distributions: (a) (u0 u0 )1/2 along the centerline and (b) (w0 w0 )1/2 along the centerline.

For visualization and identification of instantaneous vortical structures, we apply the Q criterion of Hunt et al. (1998). The Q criterion defines the regions where the second invariant of velocity gradient tensor Q is positive as vortex tubes. The quantity Q is defined by Q ¼ ðOij Oij  S ij S ij Þ=240

Fig. 11. Distribution of averaged pressure coefficient along the surface of the square cylinder.

necessary for averaging. Jeong (2003) pointed out that the required computation time necessary to get converged pressure coefficient was longer as the control parameter, fk, became smaller. He also commented that the time for at least 30 vortex shedding cycles was needed to get the converged field when fk was 0.4. In our PANS simulation, we also confirmed this tendency. Fig. 12(a) shows the time averaged pressure coefficient (CP) around the corner point B over the time covering about 25 vortex shedding cycles in the present simulation. We note the curve oscillated significantly for the case of A1. With longer computation time, this undesirable oscillation is seen to disappear to result in a smooth Cp curve as shown in Fig. 12(b). An important characteristic of the square cylinder flow is the vortex shedding. To examine how this large scale unsteadiness is predicted by the present simulations, we carry out spectral analysis of vertical velocity time series at several downstream stations. Fig. 13 shows the dominant St (fD/UN) of the present PANS simulations at x/D ¼ 3, z/D ¼ 0. The Strouhal numbers obtained in the present work are very close to the experimental data as given in Table 3. We note here again that the dominant Strouhal number is not sensitive to the grid system as can be seen in Fig. 13 and Table 3.

(18)

where Sij ¼ (ui,j+uj,i)/2 and Oij ¼ (ui,juj,i)/2. Fig. 14 shows the instantaneous coherent structures with the iso-surface of Q ¼ 0.05. All of the PANS simulations generated similar large scale coherent structures in accordance with the dominant Strouhal number predictions. As the grid system becomes finer, more realistic structures with smaller scale structures are captured as expected. Fig. 15 presents the fk distribution contour taken at the instant when the lift coefficient (CL) is smallest for various cases of the present simulation. It is seen that the instantaneous flow structures are well visualized by the contour plot of fk as fk reflects the magnitude of unresolved scales. As the grid resolution becomes finer, fk values become smaller and hence the fk contour illustrates more refined structures. Since the flow is unsteady, fk at a given position varies with time. Fig. 16 presents the time history of fk at two different positions (point 1: x/D ¼ 1.5, y/D ¼ 2.0, z/D ¼ 0; point 2: x/D ¼ 3, y/D ¼ 2.0, z/D ¼ 0.5). Point 1 is located in the center of the wake while point 2 is in the intermittent shear layer. It is seen that the fk value at point 1 ranges between 0.2 and 0.3 (case B) or 0.3 and 0.4 (case C) whereas fk at point 2 varies widely and periodically jumps to the value 1 indicating that the intermittent layer moves to RANS region (fk ¼ 1). We comment here that the time between this periodic jump of fk corresponds to the half of the Strouhal number.

5. Conclusion A two stage PANS simulation of the flow past a square cylinder was carried out. A preliminary RANS simulation was first performed to estimate the turbulent and the Kolmogorov length scales. These scales were used to devise the grid system and to initially determine the control parameter, fk. As PANS simulation proceeds, fk is updated by the equation proposed in this work. To examine the effect of grid system on the simulated results, five different grid systems were tested. These different grid systems were designed based on the turbulent length scale and Kolmogorov length scale distributions from the preliminary RANS

ARTICLE IN PRESS C.-S. Song, S.-O. Park / J. Wind Eng. Ind. Aerodyn. 97 (2009) 37–47

45

Fig. 12. Time averaged Cp distributions around the corner point: (a) averaged over 25 vortex shedding cycles and (b) averaged over 35 vortex shedding cycles.

Fig. 13. Dominant Strouhal number at x/D ¼ 3 and z/D ¼ 0: (a) Case A1, (b) Case A2, (c) Case B, and (d) Case C.

simulation. For this particular flow, it was confirmed that the Strouhal number and the drag coefficient were relatively insensitive to the grid system whereas the recirculation length was most sensitive. Comparisons of the predicted results of the present PANS simulation indicated that the grid system which generates the control parameter to be less than about 0.5 within the fully turbulent wake region produce acceptable solutions for practical purpose. This can be an effective guide to determine the spanwise grid spacing as exemplified in various simulation cases. In

turbulent flow simulations, the spanwise grid spacing is usually greater than the grid spacing in the other two directions. Thus the cut-off grid length (D) is largely determined by the spanwise grid spacing. The grid spacing for fk ¼ 0.5 corresponds to D ¼ 40Zavg (Fig. 7) which requires about 34 mesh points over the spanwise length of 4D. An overall comparison of various data from the present PANS simulations with the other available data from various sources indicates clearly that the PANS simulation performs very well

ARTICLE IN PRESS 46

C.-S. Song, S.-O. Park / J. Wind Eng. Ind. Aerodyn. 97 (2009) 37–47

Fig. 14. Q criteria based coherent structures: (a) Case A1, (b) Case A2, (c) Case B, and (d) Case C.

Fig. 15. Distribution of instantaneous control parameter: (a) Case A1, (b) Case A2, (c) Case B, and (d) Case C.

ARTICLE IN PRESS C.-S. Song, S.-O. Park / J. Wind Eng. Ind. Aerodyn. 97 (2009) 37–47

47

Fig. 16. Time history of instantaneous control parameter: (a) Case B and (b) Case C.

even with coarser grid system. This observation together with its easier implementation into an existing RANS code substantiates that the PANS approach is a very useful alternative for practical flow computations.

Acknowledgments The computation of the present work was done by modifying a baseline code provided by Prof. George Constantinescu of University of IOWA, for which we are grateful. This work was supported by Grant no. (R01-2006-000-10520-0) from the Basic Research Program of the Korea Science and Engineering Foundation. References Abdol-Hamid, K.S., Elmiligui, A., 2005. Numerical study of high temperature jet flow using RANS/LES and PANS formulations. AIAA paper 2005-5092. Abdol-Hamid, K.S., Girimaji, S.S., 2004. A two-stage procedure toward the efficient implementation of PANS and other hybrid turbulence models. NASA TM 2004213260. Baldina, J.E., Huang, P.G., Coakley, T.J., 1997. Turbulence modeling validation, testing, and development. NASA TM-110446. Barone, M.F., Roy, C.J., 2006. Evaluation of detached eddy simulation for turbulent wake applications. AIAA Journal 44 (12), 3062–3071. Bearman, P.W., Obasaju, E.D., 1982. An experimental study of pressure fluctuations on fixed and oscillating square-section cylinders. Journal of Fluid Mechanics 119, 297–321. Bosch, G., Rodi, W., 1998. Simulation of vortex shedding past a square cylinder with different turbulence models. International Journal for Numerical Methods in Fluids 28, 601–616. Constantinescu, G., Squires, K., 2004. Numerical investigations of flow over a sphere in the subcritical and supercritical regimes. Physics of fluids 16 (5), 1449–1466. Elmiligui, A., Abdol-Hamid, K.S., Massey, S.J., Pao, S.P., 2004. Numerical study of flow past a circular cylinder using RANS, hybrid RANS/LES and PANS formulation. AIAA paper 2004-4959. Fro¨hlich, J., Terzi, D.V., 2008. Hybrid LES/RANS methods for the simulations of turbulent flows. Progress in Aerospace Science 43, 349–377. Girimaji, S.S., 2006. Partially-averaged Navier–Stokes model for turbulence: a Reynolds-averaged Navier–Stokes to direct numerical simulation bridging method. Journal of Applied Mechanics 73, 413–421. Girimaji, S., Abdol-Hamid, K.S., 2005. Partially averaged Navier Stokes model for turbulence: Implementation and Validation. AIAA paper 2005-502.

Girimaji, S., Sreenivasan, R., Jeong, E., 2003. PANS turbulence model for seamless transition between RANS, LES: fixed-point analysis and preliminary results, In: Proceedings of ASME FEDSM’03 2003 fourth ASME-JSME Joint Fluids Engineering Conferences, Honolulu, Hawai. Girimaji, S.S., Jeong, E., Srinivasan, R., 2006. Partially-averaged Navier–Stokes model for turbulence: fixed point analysis and comparison with unsteady partially averaged Navier–Stokes. Journal of Applied Mechanics 73, 422–429. Hunt, J.C.R., Wray, A.A., Moin, P., 1998. Eddies stream and convergence zones in turbulent flows. In: Proceedings of the 1998 Summer Program (CTR), pp. 193–208. Jeong, E.H., 2003. Selected problems in turbulence theory and modeling, Ph.D. Thesis, Aerospace Engineering, Texas A&M University. Kogaki, T., Kobayahi, T., Taniguchi, N., 1997. LES of flow around a square cylinder. In: Direct and Large-Eddy Simulation II, ERCOFTAC series, vol. 5. Kluwer Academic Publishers, Dordrecht, pp. 401–407. Lakshmipathy, S., Girimaji, S.S., 2006. Partially-averaged Navier–Stokes method for turbulent flow: k–o model implementation. AIAA paper 2006-119. Lee, S.S., 1997. Unsteady aerodynamics force prediction on a square cylinder using k–e turbulence model. Journal of Wind Engineering and Industrial Aerodynamics 67–68, 79–90. Lubcke, H., Schmidt, St., Rung, T., Thiele, F., 2001. Comparison of LES and RANS on bluff-body flows. Journal of Wind Engineering and Industrial Aerodynamics 89, 1471–1485. Lyn, D.A., Einav, S., Rodi, W., Park, J.H., 1995. A laser-Doppler velocimetry study of ensemble-averaged characteristics of the turbulent near wake of a square cylinder. Journal of Fluid Mechanics 304, 285–319. Nikitin, N.V., Nicoud, F., Wasistho, B., Squires, K.D., Spalart, P.R., 2000. An approach to wall modeling in large eddy simulation. Physics of Fluids 12, 1629–1632. Rodi, W., 1997. Comparison of LES and RANS calculation of the flow around bluff bodies. Journal of Wind Engineering and Industrial Aerodynamics 69–71, 55–75. Rodi, W., Ferziger, J.H., Breuer, M., Pourquie, M., 1997. Status of large eddy simulation: results of a workshop. Journal of Fluid Engineering 119, 248–262. Schmidt, St., Thiele, F., 2002. Comparison of numerical methods applied to the flow over wall-mounted cubes. International Journal of Heat and Fluid Flow 23, 330–339. Sohankar, A., Davidson, L., Norberg, C., 2000. Large eddy simulation of flow past a square cylinder: comparison of different subgrid scale models. Journal of Fluids Engineering 122, 39–47. Spalart, P.R., Allmaras, S.R., 1992. A one equation turbulence model for aerodynamics flows. AIAA paper 1992-0439. Spalart, P.R., Jou, W.H., Strelets, M., Allmaras, S.R., 1997. Comments on the feasibility of LES for wings and on a hybrid RANS/LES approach. Advances in DNS /LES: First AFOSR International Conference on DNS/LES, Greyden, Columbus. Speziale, C.G., 1998. Turbulence modelling for time-dependent RANS and VLES: a review. AIAA Journal 38 (2), 173–183. Voke, P.R., 1997. Flow past a square cylinder: test case LES2. In: Direct and LargeEddy Simulation II, ERCOFTAC series, vol. 5. Kluwer Academic Publishers, Dordrecht, pp. 355–373.