See discussions, stats, and author profiles for this publication at: https://www.researchgate.

net/publication/237077897

GRID CONVERGENCE STUDY FOR A TWO-DIMENSIONAL SIMULATION OF

FLOW AROUND A SQUARE CYLINDER AT A LOW REYNOLDS NUMBER

Article · November 2009

CITATIONS READS

49 2,060

3 authors:

Mohamed Sukri Mat Ali Con Doolan

Universiti Teknologi Malaysia UNSW Sydney
88 PUBLICATIONS 523 CITATIONS 266 PUBLICATIONS 2,702 CITATIONS

SEE PROFILE SEE PROFILE

Vincent Wheatley
The University of Queensland
119 PUBLICATIONS 1,152 CITATIONS

SEE PROFILE

Some of the authors of this publication are also working on these related projects:

Noise generation of wall-mounted airfoils View project

Wind Farm Noise View project

All content following this page was uploaded by Mohamed Sukri Mat Ali on 03 April 2014.

The user has requested enhancement of the downloaded file.

 Seventh International Conference on CFD in the Minerals and Process Industries
CSIRO, Melbourne, Australia
9-11 December 2009

GRID CONVERGENCE STUDY FOR A TWO-DIMENSIONAL SIMULATION OF

FLOW AROUND A SQUARE CYLINDER AT A LOW REYNOLDS NUMBER
Mohamed Sukri Mat ALI∗ , Con J. DOOLAN∗ and Vincent WHEATLEY♠

∗
School of Mechanical Engineering, The University of Adelaide, South Australia 5005, AUSTRALIA
♠
School of Engineering, The University of Queensland, Queensland 4072, AUSTRALIA

ABSTRACT The DNS results from the three examples show

This paper describes a systematic method of refining a similar physical flow features, such as vortex shedding,
computational grid for the Direct Numerical Simulation but there are small discrepancies in the integrated
(DNS) of flow about a square cylinder. The grid variables e.g., Strouhal number (St), root mean square lift
refinement method involves two stages. The first stage coefficient (CLrms) and mean drag coefficient (CDmean). It
constructs computational meshes based on reasonable is well accepted by the CFD community that the error
estimates of cell size and grid stretching ratios and from the numerical simulation is not solely from the grid
investigates the sensitivity of the computed flow field to convergence error, but also from many other error sources.
these parameters. The findings from this stage are used as However, one can minimise the total error by reducing the
a guide for further grid refinement. In the second stage, error due to grid dependence and this must be done in a
the grid is divided into four regions to minimise the grid systematic manner.
stretching ratio. The smallest cell is sized so the boundary The purpose of this study is to assess the DNS grid
layer on the front face of the cylinder is adequately convergence of flow around a square cylinder. While the
resolved. Solutions on seven different grids are presented results of grid refinement studies are generally reported,
to investigate the effect of numerical scheme, boundary there is a need for a detailed assessment to be published to
conditions and grid independence. The level of grid assist those involved in numerical simulation. Further, for
independence is evaluated using a form of Richardson the particular case of the square cylinder, there needs to be
extrapolation and the study shows that the finest grid an objective method of determining adequate resolution of
solution has a Grid Convergence Index (GCI) of less than the boundary layers on the upstream face. This paper
5%. satisfies both these needs.
Two stages of grid refinement are involved in this
study and the level of grid convergence for the later stage
INTRODUCTION
is evaluated using the Grid Convergence Index (GCI). A
The study of flow around a square cylinder has been and novel of assessing grid independence using Thwaites’
continues to be the subject of intense research (Sohankar method is also presented. DNS results from the grid with
et al.; 1998, Saha et al.;2002, Ozgoren;2005, Inoue et al.; the lowest GCI are then compared with previously
2006 and Doolan;2009). Owing to its simplicity, this published data.
shape has become important for the study of the This report is organized as follows. After explaining
fundamental properties of bluff body wakes. These types the test case and the solution methodology, the first grid
of wakes are almost ubiquitous in all fields of engineering. refinement study is presented and suggestions for grid
For example, many types of flow mixers rely on bluff improvements are made. A second grid refinement study
body wakes to enhance mixing. is then described and the level of grid convergence is
Within the scope of predicting the flow field around a evaluated. Finally, conclusions are given in the last
square cylinder using numerical analysis, many similar section.
investigations have been made, but the results always
show small discrepancies even though the overall global
MODEL DESCRIPTION
trends are similar. One of the reasons for these
discrepancies is the difference in the construction of the Test case
mesh. Sohankar et al. (1998) for example, constructed a
non-uniform mesh near the wall of a square cylinder and, D
at 5 diameters away from the cylinder surface, a uniform
mesh was applied. The size of the smallest cell was
located at the leading edge of the square cylinder with cell SQUARE
size of h/D = 0.004. Another similar study was carried out U∞ CYLINDER D
by Inoue et al. (2006). They constructed a non-uniform
mesh but divided the computational domain into three
regions, each with a different grid ratio. The smallest cell
was located along the edges of the square cylinder with
the value of h/D = 0.01. The most recent study by Doolan Figure 1: Schematic of flow geometry.
(2009) investigated the grid convergence for three
different grid resolutions on DNS around a square cylinder
and found that the solution converged when the smallest
cell size along the square cylinder edge was h/D = 0.0167.

Copyright © 2009 CSIRO Australia 1

 CASE Nx x Ny Gedge gu gd gtb Δu Δd Δtb
A 625 x 230 100 x 4edges 1.0212 1.0072 1.0212 0.0592x0.01 0.01x0.045 0.01x0.0592
B 550 x 200 200 x 4edges 1.1142 1.0476 1.0985 0.0103x0.01 0.01x0.009 0.005x0.0181

Table 1: Grid parameter for case A and case B. Subscript u, d and tb represent upstream, downstream and top and bottom
location respectively. The total number of nodes is Nx x Ny, number of nodes along the cylinder front edge is Gedge, g is cell
size ratio between the adjacent cell and the cell area is Δ.

The case under investigation is a rigid square cylinder ratio. The case with the coarser grid is called case A and
immersed in a two-dimensional uniform incompressible the finer one is called case B (Table 1).
flow at a constant free stream velocity. Figure 1 shows a Small cells are distributed near the square cylinder
schematic illustrating the square cylinder with a side wall. The grid is then stretched from the first cell located
length of D, immersed in a flow with a constant free at the square cylinder wall to the last cell located at along
stream velocity of U∞. The Reynolds number in this study the boundary of the computational domain by adopting the
is 150 (Re = U∞ D/ν), based on the square cylinder height following equation;
D, kinematic viscosity ν and free stream velocity U∞. All
geometrical parameters and velocity flow fields are g = G 1 / ( N −1) (1)
normalized by the square cylinder height D and free
stream velocity U∞ respectively. where g is cell size ratio between adjacent cells, G is the
Solution methodology
cell size ratio between the largest and the smallest cell and
N is number of cells. Exception is made for the grid
Numerical investigation by Sohankar et al. (1999) has distribution behind the square cylinder as the wake
shown that for Re = 150 flow around a single square formation length is located at around one diameter
cylinder, two- and three-dimensional simulations give downstream from the trailing edge of the cylinder
identical results. Additionally, experimental (Luo et al.; (Zdravkovich; 1987). Therefore, the grid is made uniform
2003 and Luo et al.; 2007) and numerical (Saha et al.; until two diameters downstream from the downstream
2003) studies have shown that three dimensional edge of the cylinder before it is also stretched towards the
instability (mode A) first occurs at Re = 160. Therefore, to boundary outlet. Table 1 summarizes the grid parameters
reduce computational cost, the flow is mathematically for case A and case B and Figure 2 shows the cell
modelled in two-dimensional form. distribution around the square cylinder for case A (top)
The primitive variables of the flow fields are and case B (bottom) respectively.
calculated numerically based on the two-dimensional
unsteady incompressible Navier-Stokes and continuity Flow visualisation
equations. The OpenFOAM (Weller et al.; 1998)
numerical simulation system is used to solve these In the first stage of the grid refinement study, flow
governing equations. visualisations of each case were compared. The vorticity
The pressure implicit split operator (PISO) solution contour should shows a periodic staggered arrangement
algorithm (Barton; 1998) with two correction steps for vortices trail behind the cylinder which is known as a von
pressure-velocity coupling are used to solve these Kármán vortex street.
transient problems. The convergence criterion for pressure Figure 3 shows the instantaneous spanwise vorticity
and velocity solutions are set so the residual falls below contours for case A (top) and case B (bottom). In both
the tolerance of 10−6 and 10−5 respectively at each time cases, similar vortex shedding is observed, where the free
step. A rectilinear grid system is applied for all cases. shear layers from the top and bottom sides of the cylinder
In the first stage of the grid refinement study, the 1st- roll-up downstream of the cylinder. The roll-up continues
order implicit Euler method is used for temporal and then staggered eddies are shed downstream.
discretisation, the convection and viscous terms are However, there are two obvious discrepancies in the
discretised using the 2nd-order unbounded Gauss linear vorticity distributions. First, at the square cylinder leading
differencing scheme (central differencing). The time step edges, case A has a small scattered vortices. Meanwhile
for pressure, convection and diffusion terms is set at for case B, the vortex smoothly sheds away from the
ΔtU∞/D = 0.005 with the requirement to keep the CFL leading edge. The second discrepancy is in the wake at a
number below 0.5. few diameters downstream from the cylinder. For case A,
For second stage of the grid refinement study, the the wake shows a staggered arrangement of vortices with
temporal discretisation is advanced using a 2nd-order concentrated, well defined vortices of alternating sign
backward scheme. The convection term is discretised extending downstream. For case B, the well-defined
using a 2nd-order upwind scheme and for the viscous term, vortex structures are lost. Additionally, both cases have
a 2nd-order unbounded Gauss linear differencing scheme is small spots of vorticity in the wake with case B having
used. Three different time steps have been imposed more spots than case A.
corresponding to the three different grid cases so the CFL These observations can be related to the cell size.
number is always less than 0.5. Case A has a coarser grid around the leading sharp edges
of the square cylinder when compared with the grid of
case B. From the previous study of Sohankar et al. (1997)
GRID REFINEMENT STAGE 1
the leading sharp edges are the points where the boundary
The first stage of the grid refinement is investigated by layer starts to separate and the flow field gradients are
comparing the results from two DNS simulations whose high at these points. The vortex spots generated at the
only differences are the cell size and the grid stretching upstream sharp edges of the case A are believed due to the

Copyright © 2009 CSIRO Australia 2

large cell size. Then the simulation is unable to accurately GRID REFINEMENT STAGE 2
resolve the high velocity gradient at the leading edges. In the second stage of the grid refinement study, three
A similar reason can be used for the loss of well- different grid resolutions are used. Case C has the coarsest
defined vortex structures in case B. The high grading ratio resolution, case D has medium resolution and case E has
(gtb= 1.0985) used to stretch the grid in a transverse the finest grid. Each of the grids has four zones. Extending
direction in case B created a large discrepancy in size one diameter from the cylinder wall is region (1). In this
between adjacent cells. Thus, case B had large cells only a region, the cell size is uniform and is the overall
few diameters away from the cylinder. This type of mesh minimum. From one to three diameters from the cylinder
produces high numerical dissipation when compared with wall is region (2). The cell distribution in this region is
the cell grid size used in case A. non-uniform with the cells being stretched slowly away
from the cylinder. Region (3) extends from three to ten
diameters from the cylinder wall. The cells in this region
A) Case A were stretched and the cell size is larger than the previous
two regions. From ten diameters downstream of the
cylinder to the outlet of the computational domain is
region (4). This region is constructed to minimize the grid
stretching ratio along the wake.
The upstream, top and bottom of the boundary
computational domain are 10 diameters away from the
body and the outlet is 20 diameters downstream from the
body. This computational domain size is similar as in
simulation works by Doolan (2009) and larger than in
simulation works by Sohankar et al. (1999).
Boundary layer thickness prediction
B) Case B When fluid flows past a bluff body, boundary layers
develop from each side of the stagnation point where their
characteristics length scale can be very small. The
characteristics of the boundary layer play an important
role to the overall flow structure (i.e. flow separation),
therefore it is necessary to accurately capture the
boundary layer profile. This can be achieved by first
evaluating the boundary layer thickness based on
Thwaites' method (Thwaites; 1949).
An empirical relation for the momentum thickness (θ)
based on Thwaites integral solution for incompressible
laminar flow with a stagnation point on the body is;
Figure 2: Cell distribution around the square
cylinder. 2 L
⎛θ⎞ 0.45υ

A) Case A
⎜ ⎟ RL = 6
L
⎝ ⎠ ue
0
∫
u e5 ds (2)

where RL is defined as U ∞ L / υ and υ is kinematic

viscosity. The edge velocity, ue is calculated numerically
using potential flow solution and using OpenFOAM’s
potential flow solver to solve the equation. The integration
in eqn. (2) is carried out from the stagnation point (y=0) to
the front sharp corner of the square cylinder (y=L). Then,
the momentum thickness is transformed into the
displacement thickness using the following relationship;
B) Case B

δ * = θH(λ) (3)

where H(λ) is a shape factor. (Refer to Cebeci and

Bradshaw (1977) for the complete equation.) For laminar
flow, the following approximation for the boundary layer
thickness can be made;

δ ≈ 3δ * (4)
Figure 3: Instantaneous streamwise vorticity contours
colored by rotation direction (blue: clockwise and green:
anticlockwise), 40 equally spaced over the range - Figure 4 shows the comparison of boundary layer
4≤ΩD/U∞≤4 (Ω is spanwise vorticity). thickness (δ/D) from Thwaites' method and simulation
results of case A and case B. All three lines exhibited a
similar boundary layer thickness distribution. For the

Copyright © 2009 CSIRO Australia 3

simulation results (case A and case B), the boundary layer The influence of the outlet boundary condition
thickness is larger than the calculated value using A comparative study is carried out between two common
Thwaites's method. types of computational outlet conditions. The
conventional outlet boundary condition is defined as a
fixed pressure value and zero streamwise velocity
gradient. The second type of outlet is known as convective
boundary condition. This was proposed by Orlanski
(1976) for a problem of flow with a wavy structure. The
following equation is used:

∂U i ∂U i
+U c =0 (5)
∂t ∂xi

where indices i = 1, 2 refer to the streamwise (x) and

crosswise (y) directions, respectively and Uc is the
Figure 4: Boundary layer thicknesses from Thwaites' convective velocity. Following (Le et al.; 1997), the value
method and cases A and B. is set to Uc = 0.8U∞.
A DNS similar to case C is carried out and it is labelled
Ideally, the smallest cell size should be considerably as case F. For case C, constant pressure and zero velocity
less than the minimum boundary layer thickness. gradient are applied at the outlet, meanwhile for case F, a
However, as the smallest thickness of the boundary layer convective boundary condition is introduced. Other
from Thwaites’ method is nearly zero, it is difficult to parameters are kept unchanged. Both cases exhibited
achieve this. Additionally, potential flow theory does not similar vortex shedding patterns and the vortices pass the
able to accurately predict the rotational flow at the sharp outlet boundary smoothly. Inspection of the global results
corner. The true boundary layer thickness in this region is (CLrms, CDmean, and St) shows that both type of outlet
much larger as indicated by the simulations. Therefore, to boundary condition produced identical results.
make a practical grid, the cell size is based on the
The influence of temporal discretisation schemes
boundary layer thickness at position 10% before the sharp
corner. From Thwaites' method calculation, the boundary Two temporal discretisation schemes are compared to
layer thickness is approximately δ/D = 0.045 at that investigate the influence of higher order accuracy on the
location, y/D = 0.45 from the stagnation point. DNS results. Case D uses a 2nd-order backward Euler
Using this approximate minimum boundary layer scheme for temporal discretisation. A second case labelled
thickness, three different grids were constructed. The as case G uses a 1st-order Euler implicit method in time.
coarsest grid has approximately one cell within the Other parameters are keep identical between the two
boundary layer at that point, the medium grid has three cases.
and the finest one has five cells. Table 2 summarizes
parameters used for case C, D and E. For each case, the CASE CLrms CDmean St
grid stretching ratio in the successive regions is set as
case D 0.2917 1.484 0.1577
g (1) ≤ g (2 ) ≤ g (3) so the cells are only allowed to stretch
case G 0.2849 1.478 0.1556
slowly away from the cylinder.
Table 3: Comparison of global results between DNS with
CASE C D E 2nd-order accurate (case D) and 1st-order accurate (case G)
Nx x Ny 259 x 252 400 x 300 520 x 440 temporal discretisation.

Gedge 36 x 4edges 60 x 4edges 100 x 4edges Inspections of the global results in Table 3 shows that
the higher order scheme does not heavily influence the
Δu = Δd = Δtb 0.0278x0.0278 0.0167x0.0167 0.01x0.01 DNS results. Using a 1st-order scheme only slightly under
g(1) 1 1 1 predicts the global results when compared with the 2nd-
order temporal scheme.
g(2) 1.047 1.033 1.042
Grid convergence study using Richardson
g(3) 1.047 1.057 1.057
extrapolation
∆tU∞/D 0.002 0.005 0.005 Richardson extrapolation (Richardson et al.: 1927) is used
to calculate a higher-order estimate of the flow fields from
Table 2: Grid parameter for case C, D and case E. a series of lower-order discrete values (f1, f2, …, fn). For
Subscript u, d and tb represent upstream, downstream and the case of grid refinement study, the value estimated
top and bottom location respectively and subcript (1), (2) from the Richardson extrapolation is the value that would
and (3) represent region 1, 2 and 3 of computational results if the cell grid size tended to zero, (h→0). The
domain respectively. The total number of nodes is Nx x Ny, extrapolation is made from the results of at least two
number of nodes along the cylinder front edge is Gedge, g different grid solutions. A convergence study requires a
is cell size ratio between the adjacent cell and the cell area minimum of three grid solutions (Stern et al.: 2001).
is Δ. Roache (1994) generalized Richardson extrapolation
by introducing the pth-order methods;

[ ]
f exact ≈ f1 + (f1 − f 2 ) /(r p − 1 ) (6)

Copyright © 2009 CSIRO Australia 4

In this study, the grid refinement ratio r, is constant;

r = ΔC / ΔD = ΔD / ΔE = 1.667 (7)

From equation (6), the extrapolated value is varied by

different choice the order p. According to Stern et al.
(2001) the order-of-accuracy can be estimated by using
the following equation;
ln(ε 32 / ε 21 )
p= (8)
ln(r)

ε i+1,i = f i+1 − f i (9)

Figure 5: Comparison of global results normalized by

To evaluate the extrapolated value from these solutions,
the extrapolated value, between three grid solutions and
the convergence conditions of the system must be first
Richardson Extrapolation estimation.
determined. The possible convergence conditions are;
CLrms; --□-- and CDmean; --◊--
1. Monotonic convergence; 0<R<1
As listed in Table 4, there is a reduction GCI value
2. Oscillatory convergence; R<0
for the successive grid refinements (GCI21 < GCI32) about
3. Divergence; R>1
each of the three variables. The GCI for finer grid (GCI21)
is relatively low if compared to the coarser grid (GCI32),
where R is the convergence ratio and it is determined by
indicating that the dependency of the numerical simulation
the equation;
on the cell size has been reduced. Additionally, as the GCI
ε
R = 21 (10) reduction from the coarser grid to the finer grid is
ε32 relatively high, the grid independent solution can be said
to have been nearly achieved. Further refinement of the
Table 4 summarizes the order of accuracy for root grid will not give much change in the simulation results.
mean square lift coefficient (CLrms) and mean drag For variable CLrms, the extrapolated value is only
coefficient (CDmean) from the simulation results on three slightly lower than the finer grid solution (h/D = 0.01) and
grids. All means and root mean square results were it is in the range of the finer grid GCI as shown in figure
obtained once a statistically stationary flow was 5. Similar behavior is observed for variable of mean drag
observered. The convergence condition for CLrms and coefficient. Therefore, it is shown that the solution has
CDmean are monotonic. Strouhal number (St) is not converged with the refinement from the coarser grid to the
included in this analysis as the differences between the finer grid.
successive grid solutions are very small (St = 0.1556, Similar behavior is observed for the point variable of
0.1577 and 0.1602 for case C, D and E respectively). boundary layer thickness (δ/D). The boundary layer
Hence, this parameter is not suitable for use in this grid thickness is calculated at the leading edge of the square
refinement study. cylinder. As the grid is refined, the discrepancy between
The Grid Convergence Index (GCI) provides a the solution and the extrapolated value become small. In
uniform measure of convergence for grid refinement this paper, the discrepancy between the simulation value
studies (Roache:1994). It is based on estimated fractional and this extrapolated value is used to define the error,
error derived from the generalization of Richardson f − f RE
extrapolation. The GCI value represents the resolution Ei = i (12)
f RE
level and how much the solution approaches the
Table 5 shows that the convergence condition is
asymptotic value. The GCI for the fine grid solution can
monotonic. The order-of-accuracy is then calculated by
be written as;
using equation (8) and (9) using the L2-Norm of the
GCI i+1,i = Fs
|εi+1,i | (11) values. Figure 6 shows that the successive grid
f i (r p − 1 ) refinements nearly achieved the asymptotic value at the
Following Wilcox (2006), the safety factor (Fs) selected finest grid resolution, where the relative error compared
for this study is 1.25. with the RE is only E1 = 0.28%. This indicates that the
boundary layer has been adequately resolved.

ε32 ε21 GCI32 GCI21 Variable ||ε32/f1|| ||ε21/f2|| <R> p GCI32 GCI21
R p
(10-2) (10-2) (%) (%)
δ 0.0050 0.0172 0.275 2,(2.54) 1.21% 0.35%
CLrms 1.53 0.64 0.4183 1.7062 4.71 2.02
Table 5: Order of accuracy and Grid Convergence Index
CDmean 2.30 1.0 0.4348 1.6305 1.49 0.65
for point variable. Subscripts 3, 2 and 1 represent case C,
D and E respectively. The order of accuracy p for both
Table 4: Order of accuracy and Grid Convergence Index variables are taken as 2 as the calculated values (indicated
for three integration variables. Subscripts 3, 2 and 1 in a brace) exceeded the order of accuracy of the
represent case C, D and E, respectively. discretisation scheme.

Copyright © 2009 CSIRO Australia 5

 further DNS analysis as the GCIs for all variables being
studied were less than 5%.

REFERENCES
BARTON, I. E. (1998). “ Comparison of SIMPLE- and
PISO- type algorithms for transient flows.” International Journal
for Numerical Methods in Fluids, 26,459-483.
CEBECI T. and BRADSHAW P. (1977). Momentum
Transfer in Boundary Layers. McGraw-Hill Inc.,US.
DOOLAN, C. J. (2009). “Flat-plate interaction with the near
wake of a square cylinder.” AIAA Journal, 47,475-478.
INOUE, O., IWAKAMI, W., and HATAKEYAMA, N.
(2006). “Aeolian tones radiated from flow past two square
cylinders in a side-by-side arrangement.” Physics of Fluids,
18(4),046104.
Figure 6: Percent of RMS error for boundary layer LE, H., MOIN, P., and KIM, J. (1997). “Direct numerical
thickness along the leading edge. The order-of- accuracy, simulation of turbulent flow over a backward-facing step.”
p = 2 and E1 = 0.28%, E2 = 0.77% and E3 = 1.70%. Journal of Fluid Mechanics, 330,349-374.
LUO, S. C., CHEW, Y. T., and NG, Y. T. (2003).
The grid convergence index (GCI) is estimated by using “Characteristics of square cylinder wake transition flows.”
equation (11) with a safety factor of Fs = 1.25. As listed in Physics of Fluids, 15(9),2549-2559.
the same table (Table 5), successive grid refinements LUO, S. C., TONG, X. H., and KHOO, B. C. (2007).
resulted in a GCI reduction. Therefore, that can be said “Transition phenomena in the wake of a square cylinder.” Journal
that the solution on the finest grid resolution is nearly grid of Fluids and Structures, 23(2),227-248.
independent. ORLANSKI, I. (1976). “A simple boundary condition for
unbounded hyperbolic flows.” Computational Physics, 21(3),251-
Comparison of numerical data from case E with 269.
previous studies OZGOREN M. (2005). “ Flow structure in the downstream of
a square and circular cylinders.” Flow Measurement and
Table 6 compares the current DNS results of case E with Instrumentation, 17,225-235
the summarized previous available data taken from RICHARDSON, L. F. and GAUNT, J. A. (1927). “The
Doolan (2009). The Strouhal number and root mean deferred approach to the limit. Part I. Single lattice. Part II.
square lift coefficient from case E are in excellent Interpenetrating lattices.” Philosophical Transactions of the Royal
agreement and in the range of the previous available data. Society of London. Series A, Containing Papers of a
Mean drag coefficient is slightly higher than previous Mathematical or Physical Character, 226,299-361.
ROACHE, P. J. (1994). “Perspective: A method for uniform
studies, but is still in good agreement.
reporting of grid refinement studies.” Journal of Fluids
Engineering, 116(3),405-41
CLrms CDmean St SAHA, A. K., BISWAS, G. and MURALIDHAR, K. (2002).
Experiments 0.148 “ Three-dimensional study of flow past a square cylinder at low
(Okijama:1982, - 1.40 to Reynolds numbers.” International Journal of Heat and Fluid
Sohankar et al.:1999) 0.155 Flow, 24,54-66.
SOHANKAR, A., NORBERG, C., and DAVIDSON, L.
DNS(Doolan;2009) 0.296 1.44 0.156 (1997). “Numerical simulation of unsteady low-reynolds number
flow around rectangular cylinders at incidence.” Journal of Wind
DNS (Sohankar et al.;1998) 0.230 1.44 0.165
Engineering and Industrial Aerodynamics, 69-71,189-201.
DNS (Inoue et al.:2006) 0.40(peak) 1.40 0.151 SOHANKAR, A., NORBERG, C., and DAVIDSON, L.
(1998). “Low-Reynolds-number flow around a square cylinder at
current DNS (case E) 0.285 1.47 0.160 incidence: study of blockage, onset of vortex shedding and outlet
boundary condition.” International Journal of Numerical
Table 6: Comparison of current DNS (case E) with previous Methods in Fluids, 26(1),39-56.
studies listed by Doolan (2009). SOHANKAR, A., NORBERG, C., and DAVIDSON, L.
(1999). “Simulation of three-dimensional flow around a square
III. CONCLUSION cylinder at moderate reynolds numbers.” Physics of Fluids,
11(2),288-306.
For the first grid refinement study, flow visualization of STERN, F., WILSON, R. V., COLEMAN, H. W., and
vorticity contours from the two different grid resolutions Paterson, E. G. (2001). “Comprehensive approach to verication
exhibited similar behavior of some vortex spots a few and validation of cfd simulations part 1: Methodology and
diameters downstream from the square cylinder. procedures.” Journal of Fluids Engineering, 123(4),793-802.
THWAITES, B. (1949) “Approximate calculation of the
Additionally, for case A there was a small scattered vortex
laminar boundary layer” Aeronautical Quarterly, 1,245-280.
near the upstream corner of the square cylinder. These WELLER, H. G., TABOR, G., JASAK, H., and Fureby, C.
observations indicated that it was necessary to refine the (1998). “A tensorial approach to computational continum
grid in these regions and the grid stretching ratio between mechanics using object-oriented techniques.” Computers in
cells should be minimized. Physics, 12,620-631.
In the second grid refinement study, inspection of WILCOX, D.C. (2006), Turbulence Modeling for CFD, 3rd
GCI values for integration and points variables shows that Ed., DCW Industries, Inc..
there was a gradual reduction when the grid system was ZDRAVKOVICH, M. M. (1987). “The effects of interference
between circular cylinders in cross flow.” Journal of Fluids and
refined. Comparison between the extrapolated value
Structures, 1(2),239-261.
calculated from Richardson extrapolation indicated that
the finer grid (case E) was appropriate to be used in

Copyright © 2009 CSIRO Australia 6

View publication stats