A Numerical Investigation of Turbulent Flow Over Single and Tandem Square Cylinders

This paper presents the three-dimensional flow investigations over the square cylinders placed in the tandem arrangement. Two different flow configurations were investigated in detail; one comprising of a single square cylinder and the other comprising of three square cylinders placed in a tandem arrangement with the spacing of six times the width (w) of each square cylinder. The Reynolds number based on the width of the square cylinder (w) and free stream velocity (U o ) is 22,000. The problem was solved numerically using an Unsteady Reynolds-Averaged Navier-Stokes (URANS) based model and Large Eddy Simulation (LES) based model. Strouhal Number, lift, and drag coefficients were computed for each configuration. By comparing both the models using contour plots of pressure, velocity and vorticity it is found that the LES model is more accurate to capture the turbulence around single and tandem square cylinders. In the tandem arrangement, complex periodic vortex shedding was observed in the wake of each square cylinder. The production of turbulent kinetic energy was also investigated to understand the roles of stresses during flow over the cylinders. The analysis showed that the production of turbulence by normal stresses is higher than that of shear stresses. Furthermore, it was observed that the flow over the first cylinder arranged in tandem is quite identical to that of the single square cylinder. Moreover, the upstream cylinder experienced a higher lift in comparison to the downstream cylinders.


INTRODUCTION
ver the last century, flow over bluff bodies is the subject of intense investigations both numerically and experimentally. The flow over bluff bodies is known for complexities including flow separation, wake generation, vortex shedding, and flow-induced vibrations. These flows have importance in electrical transmission lines, bluff shaped tall buildings like skyscrapers in a city, stacks, research is highlighted in this work.
Williamson [1] analyzed the vortex dynamics in the wake of a circular cylinder over a wide range of Reynolds numbers (Re). It was investigated that the natural transition to turbulence in the cylinder wake was due to the amplified vortex dislocations that are present in the near wake of the cylinder. Mode A was characterized by the appearance of stream-wise vortices with a wavelength of approximately fourcylinder diameters (4d) in the span-wise direction, and it was identified in the secondary instability of the wake. Randomly occurring dislocations seemed to appear along the span of the cylinder which triggered the turbulence. Williamson [2] also studied threedimensional (3-D) wake patterns in the laminar regime and observed fine-scale 3-D span-wise instabilities and the appearance of vortex dislocations in the wake transition regime. The new oblique wave resonance mechanism in the far wake region was discovered.
The Digital Particle Image Velocimetry measurements were used by Govardhan and Williamson [3] to study the mean and unstable velocity fields in the wake of a freely oscillating circular cylinder, at Re = 3000 to 4000. For comparison, a stationary cylinder case at approximately similar Reynolds number (Re = 3900) was also studied. It was evaluated that Reynolds stresses were significantly larger at higher Reynolds numbers but the periodic component of stress was quite similar to the stationary case. Stocks et al. [4] experimentally studied the flow field in the wake of a finite-length circular cylinder at Re of 180 to 190. It was found that natural transition to turbulence was due to an amplification of vortex dislocations that are generated in the near wake of the cylinder (mode A). They also observed that the wake properties are strongly affected by a transition in the central, downstream plane of the cylinder. They finally stated that naturally occurring turbulence can only be observed after the formation of mode A.
A detailed Floquet analysis of the transition around a circular cylinder was analyzed by Thompson et al. [5]. In this study, the growth of an elliptic instability in the forming vortex cores was observed. The amplification of a strong strain field was observed in the hyperbolic region that was in between the forming and shed vortices. It was analyzed that the wake immediately behind the cylinder showed distinct signs of a cooperative elliptic instability (counter-rotating vortices). Further downstream, after the vortices had been shed into the wake, the instability again was grown in the cores. It was concluded that the initial instability causing the transition gave rise to two shedding modes referred to as modes A and B, which led to the rapid evolution to a fully turbulent flow. Robichaux et al. [6] investigated the threedimensional instabilities of the square cylinder and showed that there were at least two instabilities that were similar to the circular cylinder. In particular, mode A and mode B three-dimensional instabilities were observed, both for the square and circular cylinders. In that study, a third mode referred to as mode S was observed in the square cylinder, which had not been previously reported. This mode was absent in the case of a circular cylinder, at least up to Re=300. The wavelength of mode S was approximated as twice of mode B and most probably observed in the rectangular-shaped bodies named as 2T-periodic three-dimensional Floquet instability.
The division of flow regimes in a square cylinder wake at various angles of attack (α) was examined by Luo et al. [7]. The wavelengths of the stream-wise vortices in the modes A and B regimes were measured by using the auto-correlation method. The critical Reynolds number (Recr) for the inception of these instability modes were identified through the determination of discontinuities in the Strouhal number (St) versus Reynolds number curves. The Particle Image Velocimetry technique was used for quantitatively measuring the parameters of wake vortices at Recr in the range from 60 to 350. Tong et al. [8] performed similar experiments using two flow facilities one of which was a vertical water tunnel and the other was a horizontal water channel with varying α.
Lyn et al. [9] experimentally studied the unsteady turbulent flow over a square cylinder at Re of 21,400 through phase-averaged velocity statistics using twocomponent Laser Doppler velocimetry measurements. Mean, periodic, and random components of the flow field over a square cylinder were obtained. Flow structures were examined and a distinction between vorticity and streamline saddles was drawn. According to this study, the vortices in near wake regions were found to be different from vortices in the far-field, where they grew to maturity. Streamline centers were identified as local maxima of vorticity and were associated with local maxima in Reynolds normal stresses (u´u´) or turbulent kinetic energy (k). The turbulent kinetic energy (k) is defined as k = (u u + v v + w w ). Streamline saddles coincide with high Reynolds shear stresses (u´v´) and high rates of turbulence production (Pk). The work of Lyn et al. [9] now stands as a benchmark problem for the numerical investigations of flow over a square cylinder at high Reynolds number.
The flow field in the wake region of different bluff bodies including circular, square and triangle crosssection cylinders was experimentally and numerically investigated by Yagmur et al. [10]. The Particle Image Velocimetry (PIV) method was used for experimental studies. Experiments were performed in an open water channel at Re of 5000 and 10000 defined according to the characteristic lengths of the cylinders. The transient simulation with LES based turbulence model from ANSYS-Fluent Software was used for the numerical work and then validated by the experiments [10]. It was examined that the numerical and experimental results had a good agreement in respect of the instantaneous and time-averaged flow field patterns of vorticity, and the velocity component in the streamwise direction.
Mahir [11] performed the 2-D and 3-D numerical simulations to investigate the effects on the flow structure around a square cylinder placed near a plane wall. The gap ratio between the cylinder and wall was varied from 0.2 to 4 for the Re of 175, 185 and 250. This analysis showed that the 3-D simulations predicted lower mean drag coefficients than those of the 2-D simulations. Moreover, the mean drag coefficients were slightly decreased at the large gap ratios, and these values were different from that of an isolated cylinder in a free stream case.
Sheard et al. [12] investigated the 3-D transition scenarios in the wake of a square cylinder with variation in the incidence angle of 0º ≤ α ≤ 45º over a range of Reynolds numbers. 3-D instability modes were predicted and computed using a linear stability analysis technique. These unstable modes were the long-spanwise-wavelength mode A and the shortspanwise-wavelength mode B and were present in the flow around the cylinder. However, a third mode, referred to as mode C, was also present depending on the position of cylinders. The transition to 3-D (turbulent) flow through either a mode A instability or a sub-harmonic mode C instability was found on varying the angle of incidence. Sheard [13] also studied the stability of the flow behind a cylinder with a square cross-section and focus on small incidence angles 0º to 12º. The first-occurring mode A instability was found to be completely suppressed as α was increased through 2º to 3º. It was examined that the switch from quasi-periodic to subharmonic properties occured as α was increased from 2º to 3º. The Reynolds numbers with varying incidence angle regimes for linear stability were comprehensively recorded.
Minguez et al. [14] presented an experimental and numerical analysis of turbulent flow past a square cylinder at Reynolds number of 21,400. For experimental investigations, the technique of Laser Doppler Velocimetry (LDV) was used. For numerical investigations, the LES method based on a spectral vanishing technique was employed. The flow was found to be separated at the leading edges of a cylinder. The occurrence of the 3-D Kelvin-Helmholtz (KH) pairings and von Kármán (VK) vortex shedding in the near wake region was also investigated based on visual and frequency analysis. Vortex-shedding by LES and compared with RANS based turbulence models was also investigated by Rodi [15]. Flow past a square cylinder was analyzed at Reynolds number of 22,000. For RANS based simulations, standard k-ε model and Reynolds-stress model were tested. It was found that the RANS based models were highly unpredictable to investigate the levels of turbulent fluctuations whereas LES could predict these fluctuations more accurately.
Recently, Zhang [16] compared different computational models including the Wray-Agarwal turbulence model (WA), Spalart-Allmaras (SA) model, the Shear Stress Transport (SST) k-w model and the standard Wilcox k-w model for circular and square cylinders. The computational Reynolds numbers were Re = 6.7×10 5 , 1×10 6 and 3.6×10 6 for circular cylinder and Re = 2.2×10 4 for a square cylinder. The WA model captured flow properties more precisely and exhibited good accuracy compared to other models (SA model and SST k-w model).
Yu and Kareem [17] numerically simulated the pressure and velocity field around the rectangular cylinder for Reynolds number of 10 5 using the LES scheme. They visualized the flow field over the rectangular cylinder by employing the streak-lines; which clearly demonstrated the shedding of two vortices. The flow over a square cylinder was also investigated using Smogorinsky (S-Model) and Dynamic One-Equation (OE-Model) LES subgrid models by Sohankar [18]. The turbulence created by the flow over a square cylinder for low Re of 1000 to high Re of 5×10 6 was studied. It was observed that at higher Re the size and strength of the small-scale structures were changed. Sohankar [19] also analyzed the 3-D unsteady flow over a sharp-edged rectangular cylinder with various side ratios (R) ranging from 0.4 to 4.0 at Re of 1×10 5 . The effect of Re ranging from 2×10 4 to 10 5 for R = 0.62 was also examined. It was found that the flow separated at the leading edges of the cylinder and the separated shear layer did not reattach itself to the sides of the cylinder. Likewise, the side surfaces were completely engulfed by the separation bubbles. Ochoa and Fueyo [20] examined flow over the square cylinder using Smagorinsky LES subgrid model at Reynolds number of 21400. A comparison between LES and RANS with the k-ε model was made to investigate the flow field and vortex shedding frequencies. It was also confirmed by this study that LES could predict the flow field more precisely than the RANS method.
Tong et al. [21] presented a numerical study for the flow over the chamfered square cylinders at varying degrees of chamfer and incident angles. It was concluded that chamfering reduced the drag coefficient compared with the sharp-cornered case. Whereas a slight decrease in drag was also observed as the Reynolds number increased. Conversely, the increase in drag was observed with the increase in the angle of attack. Recently, Olsen et al. [22] performed an intensive literature review including a study and comparison of different CFD models for analyzing the flow over the single square cylinder. They concluded that inlet turbulent intensity had less significant effects on the square sharp edge shaped bluff bodies. Furthermore, the blockage ratio of greater than 10% could significantly affect the drag coefficients.
Although there are extensive investigations of flow over a single square/ circular cylinder and tandem arrangements of circular cylinders, there are remarkably few investigations were made for the flow over a row of non-circular cylinders at high Reynolds number.
Carmo et al. [23] investigated the flow around two circular cylinders in a staggered arrangement with fixed streamwise separation of 5d and cross-stream separation varying from 0 to 3d. The Floquet stability analysis and Direct Numerical Simulations (DNS) technique were used to present the results. The wake transition was compared to that of a single isolated cylinder and initially unstable modes were captured accurately that first appeared in the wake transition of the flow around a single cylinder. These unstable modes A and B were present in the flow around the staggered arrangements. However, a third mode C was also present depending on the position of the cylinder.
Mital and Kumar [24] employed the stabilized Finite Element (FE) formulation to study flow-induced oscillations of a pair of equal-sized cylinders in tandem and staggered arrangement placed in the uniform incompressible flow. For various values of the structural frequency of the oscillator, 2-D computations were carried out for the Re of 100. The cylinders were separated by 5.5d in the streamwise direction and for the staggered arrangement the spacing between the two cylinders is 0.7d. It was found that the downstream cylinder that lied in the wake of the upstream cylinder experienced wakeinduced flutter.
Etminan et al. [25] focused on the unconfined flow characteristics in tandem square cylinders arrangement in both steady and unsteady periodic laminar flow regimes for the Re of 1 to 200. The fluctuating lift force acting on the upstream cylinder was found to be influenced by the phase of the flow pattern of the downstream cylinder. It was found that the lift force acting on the upstream cylinder maximized when the flow over the downstream cylinder was in phase with the upstream cylinder.
Rosales et al. [26] numerically investigated the unsteady flow-field and heat-transfer characteristics for a tandem pair of square cylinders in a 2-D laminar channel flow. The drag, lift, and heat transfer coefficients from the downstream heated cylinder were studied at the Re of 500 based on the side length. The results showed that the drag coefficient (CD) and the Nusselt number (Nu) decreased as the heated cylinder approached the wall. The highest Strouhal number (St) value was observed in the channel center. Complex periodic vortex shedding developed when the eddy promoter was offset from the primary cylinder. An incompressible Finite Volume computer program, Fluid and Heat Transfer Solver (FAHTSO), was used to perform the numerical calculations.
Chatterjee et al. [27] numerically investigated the flow past a row of five square cylinders placed side-by-side at different separation ratios at a Re of 150. The wake structure and the corresponding vortex shedding mechanism were studied extensively for different spacing between the cylinders. Depending on the separation ratio, the different flow patterns were observed. These flow patterns were supposed to be a consequence of the interaction between two types of frequencies, the vortex shedding (primary) and the cylinder interaction (secondary) frequencies. The lift and drag coefficient were found to be sinusoidal in nature which ensured that the primary (Strouhal) frequency pertaining to the vortex shedding was the dominant frequency and there waws no significant contribution of the secondary-frequency.
Fluid forces and their discontinuity acting on two square prisms in a tandem arrangement were experimentally analyzed by Alam et al. [28] at Re of 5.6×10 4 . In this investigation, the approaching flow was controlled by using a thin flat control plate. The maximum reduction in time-averaged drag was 84% for the upstream prism, and 66% for the downstream prism. It was observed that this reduction occured when the shear layer that was separated from the control plate reattached at the edges of the front face of the upstream prism.
Nazari et al. [29] used the k-w-υ ²-f model for three different configurations. They investigated the flow over a single square cylinder at Re of 13,000 and 22,000, the flow over cylinders placed in tandem at Re of 16,000 and 56,000, and the flow over a cube at a Re of 40,000. It was found that the flow over these bluff bodies became unsteady as per increment in the value of Reynolds number. It was concluded that the examination of coherent vortex shedding in complex flow using unsteady RANS with υ 2 -f model was reliable and cost-effective.
The numerical studies based on the Steady RANS (SRANS) approach did not capture the turbulent structure of the flow and SRANS cannot really predict any turbulent structures. Also, the predicted turbulent quantities are particularly poor [30]. Thus the Unsteady RANS (URANS) have been used to improve the performance in this turbulent flow case. However, the SRANS approach performs very poorly in terms of Reynold's stress predictions whereas the URANS approach is clearly much better and showed much improvement over the SRANS approach and produces results with a much better agreement with the experimental data [30]. However, even the URANS approach fails to predict Reynolds stresses accurately which is mainly due to the fact that URANS may be able to capture some large-scale unsteady flow motions but it cannot capture turbulence well. Moreover, when the detailed turbulent information is needed then it is only possible to employ LES [30,31].
The literature review shows that intensive research related to a single circular and square cylinder is carried out but less consideration is given to flow over a row of square cylinders at high Reynolds number (Re > 10 5 ). This paper presents the investigation of flow over three square cylinders placed in the tandem arrangement at Re of 22,000. Firstly, the flow over a single square cylinder (Case-I) was investigated using the URANS method and compared with the experimental data of Lyn et al. [9]. Secondly, the flow over three square cylinders placed in a tandem arrangement (Case-II) was analyzed using URANS and LES methods. The Shear-Stress Transport (SST) k-w model was employed for the application of the URANS method. Besides that from the LES method, the Dynamic Smagorinsky Model (DSM) was used to capture the turbulence. There was a need to investigate the production of kinetic energy (k) in order to understand the role of stresses that develop during the flow over a row of square cylinders. Coherent vortex shedding is observed in the wake of each cylinder and detailed examinations are presented in the result section.
The paper starts with a detailed literature review (section-1) followed by the complete problem description in Section 2. The details of geometry and grids used for single and tandem square cylinders for both URANS and LES methods are presented in Section 2. In Section 3, the details of modeling and simulation are highlighted. The URANS based SST kw, and LES based DSM models governing equations are illustrated in Section 4.1 and Section 4.2 respectively. Finally, in Section 5 the results are elaborated case viz. Section-5.1 for Case-I and Section-5.2 for Case-II. The conclusion of the present work is presented in Section 6.

PROBLEM DESCRIPTION
The problem under investigation comprises single and three square cylinders aligned with their axes perpendicular to incoming free streamflow. The schematic diagram of the domain used for the simulation of a single square cylinder of size w exposed to free stream velocity (Uo) is shown in Fig.   1a. The Reynolds number (Re = Uo w /ν) based on Uo and w is 22,000.
The section of the finest grid used for the flow investigation over a single square cylinder is shown in Fig. 1b. The 3-D computational domain size for the single square cylinder case is 26w × 20w × w and the number of grid points for the finest, medium and coarse grids are 323×259×49, 251×211×33 and 187×163×17, respectively. Grid sizes used for the single square cylinder (Case-I) and the minimum distances from cylinder surfaces to the nearest grid point are summarized in Table 1. The 3-D computational domain size for the single square cylinder case is based on the study of Sohankar et al. [18,19].
The computational domain size for URANS simulations for the three square cylinders Case-II is 40w × 16w × w and the numbers of grid elements in x, y, and z directions are 551×243×41 for the finest grid, 431×203×25 for the medium grid and 287×139×25 for the coarse grid. Grid sizes used for tandem square cylinders Case-II with the details of minimum distances from cylinders surfaces to the nearest grid points are summarized in Table 2.    For LES simulations of Case-II, the geometrical domain used is shown in Fig. 2a in which three square cylinders of size w are exposed to free stream velocity Uo. The distance between the square cylinders is six times the width (6w). The Re based on free stream velocity (Uo) and width (w) is 22,000. The section of the grid used for LES is shown in Fig. 2b. The 3-D computational domain used for LES investigation is 40w × 16w × 4w that is selected on the basis of a study conducted by Sohankar et al. [33]. The mesh comprises of 527×165×73 elements in x, y and z directions respectively ( Table 2). The problem discussed involves complex flow behavior and to capture that the grids used have y+< 1, thus, it is expected that turbulence at near-wall regions is accurately captured. ICEM CFD was used for 3-D geometry preparation and grid generation. FLUENT used unstructured grids while the structured grid was exported for FASTEST.

URANS
First, the flow over single cylinder Case-I was investigated. Later the flow over three square cylinders Case-II was simulated. In both cases, the commercial CFD code Fluent from ANSYS Inc. (version 12) was used. The code was based on the finite volume technique for spatial discretization.
Computations were performed on hexahedral grids using the SST k-ω model. Second-order upwind discretization along with the SIMPLE pressurevelocity coupling [34] was used for time integration, and adaptive time-stepping was applied to obtain timeaccurate simulations.
All computations were carried out at the Computational Thermo-fluid Laboratory (NED University of Engineering and Technology) on an HP-Pro-Liant Cluster using eight AMD Quad Core processors each with 8 GB RAM for simulations. Flows were evolved up to 526 seconds for Case-I, and 583 seconds for Case-II, which takes approximately 1536 CPU hours for a single cylinder and 6144 CPU hours for tandem arrangement, respectively.

LES
The turbulent flow field for Case-II was also investigated using the Dynamic Smagorinsky model (DSM-LES method) proposed by Germano et al. [35]. Germano et al. [35] proposed a procedure that allows the determination of the Smagorinsky constant dynamically. Through this procedure the Smagorinsky constant will be no longer a constant, but rather become more realistically a function of space and time. The computations were conducted using a noncommercial fully parallel CFD code, FASTEST. This code was based on a finite volume method and solves he 3-D filtered Navier-Stokes equations on blockstructured hexahedral meshes, using the cell-centred collocated variable arrangement. The coupling of the velocity and pressure fields was achieved with the SIMPLE algorithm. The convective transport of all variables was discretized by a second-order central differencing scheme. The second-order implicit Crank-Nicolson method was employed for time discretization.
For LES, sixteen AMD Quad-Core processors, each with 8 GB RAM were used for the simulations. LES computations were also performed at Computational Thermo-fluid Laboratory on HP-Pro-Liant Cluster with virtualization. The total simulation time for LES was 16896 CPU hours. Simulation results were averaged over 20 periods of vortex shedding, which was found sufficient to achieve the statistical significance of the results.

The Boundary Conditions (BC) for URANS and LES
A uniform velocity profile was specified at the inlet of the domain, as prescribed by Sohankar [19]. Whereas, the turbulence intensity (I) is calculated from I=0.16×(Re -1/8 )=4.51% and used as the inlet profile turbulence parameter [15,36]. In the case of LES, a convective boundary condition (BC) was prescribed at the outlet whereas, in URANS simulations, pressure outlet BC was chosen as the Dirichlet BC [37]. A periodic BC opted in the spanwise direction. The noslip BC was selected for cylinder walls [18] and at the top and bottom walls of the domain. Furthermore, it was observed that a boundary layer of thickness δ/w= 0.1, developed at the top and bottom wall of the domain and it did not affect the flow over the square cylinders. Blockage parameter β=w/H for the domain was selected on the bases of the study conducted by Sohankar [18]. Here H is the spanwise width of the domain. The value of β for URANS is 0.05 and for LES it is 0.0625.

Shear Stress Transport (SST) k-ω Model
The Reynolds-Averaged Navier-Stokes equations for unsteady, incompressible, Newtonian fluid given by Hinze [38,39] is used in this study. The continuity equation is given as: and the momentum equation is given as: The Reynolds-averaged approach to turbulence modelling requires that the Reynolds stresses in equation (2) are appropriately modelled. A common method employs the Boussinesq hypothesis to relate the Reynolds stresses to the mean velocity gradients: Turbulent kinetic energy (,) is defined as, In the present work, the Boussinesq hypothesis is used in the SST k-ω model. The advantage of this approach is the relatively low computational cost associated with the computation of the turbulent eddy viscosity (µt). In the case of the SST k-ω model, two additional transport equations (for the turbulence kinetic energy, k, and the specific dissipation rate, ω) is solved and µt are computed as a function of k and ω. The disadvantage of the Boussinesq hypothesis as presented is that it assumes µt is an isotropic scalar quantity, which is not strictly true. Moreover, when density variations are sufficiently small the Boussinesq approximation is valid [41].
For simulation with the SST k-ω model, the fluent solver is used. This model has the ability to effectively blend the formulation of the k-ω model in the nearwall region with the free-stream independence of the k-ε model in the far-field. These features make the SST k-ω model more accurate and reliable for a wider class of flows (e.g., flow past bluff bodies, aerofoils, etc.) than the standard k-ω model. Other modifications include the addition of a cross-diffusion term in the ω (eddies dissipation rate) equation and a blending function to ensure that the model equations behave appropriately in both the near-wall and far-field zones.
The SST k-ω model uses the following two equations [42] given below; and In these equations, G 0 -represents the generation of turbulence kinetic energy due to mean velocity gradients. G 0 3 represents the generation of ω. Γ 4 and Γ 5 represents the effective diffusivity of k and ω, respectively. Yand Y 3 represent the dissipation of k and ω due to turbulence. Whereas, the Sand S 3 are user-defined source terms; see further details in [34,40].

Large Eddy Simulations with Dynamic Smagorinsky Model (DSM)
LES is based on the premise that large eddies in the flow are dependent on the geometry while the smaller scales are more universal. Mathematically, one may think of separating the velocity field into a resolved and sub-grid part. The resolved part of the field represents the "large" eddies, while the subgrid part of the velocity represents the "small scales" whose effect on the resolved field is included through the subgridscale model. Formally, a function with a filtering kernel is used, In the above equation, x denotes space coordinates and t denotes time. The G is a normalized weight function or filter. However, most practical (and commercial) implementations of LES use the grid itself as the filter and perform no explicit filtering. The filtered equations are developed from the incompressible Navier-Stokes momentum equation. Substituting in the decomposition and then filtering the resulting equation gives the equation of motion for the resolved field, The stresses at the sub-grid level are described by the equation: Here μ @ is the sub-grid viscosity and S #$ is filtered strain rate.
In the present paper, the sub-grid model is actually a modified version of the DSM proposed by Germano [35].

Single Square Cylinder (Case-I)
In this section, the flow over a single square cylinder (named as Case-I) is discussed in detail. Table 3 summarizes the comparison of Strouhal number (St = f w/Uo) and the drag coefficient (CD =2FD/ρUo 2 A, here FD is the drag force, ρ is the mass density of the fluid, Uo is the flow speed of the object relative to the fluid, A is the reference area). It is found that the fine grid (St =0.134 and CD=2.170) shows the best results in comparison with the experimental results of Lyn et al. [9] and with the numerical results of Sohankar et al. [33], Sohankar [18], Nazari et al. [27] and Olsen et al. [22]. The Nazari et al. [29] and Rodi [15] data are used here as reference study for present results comparison only, however, they have used a 2-D grid in their study (as shown in Table 3).
For results elaboration, Fig 3(a) is used here to explain that (x'/w, y'/w = 0, 0) is the center of the square cylinder, whereas it is x/w=10.5, and y/w =10.0 away from the actual origin of the geometry. The streamlines at z/w = 0.5 for Case-I are shown in Fig.  3(b) and the data are measured and referenced from the center of the cylinder (x'/w =0, y'/w =0). Three distinct rotating streamlines regions can be observed, one at the top of the cylinder, second at the bottom and third large rotating region is developed in the near wake of the cylinder. The normalized length of recirculation zones (LR/w) for the top and bottom regions are calculated and it is found to be 0.82 and 0.89, respectively; whereas the length of recirculation zones in the wake region of the cylinder is 1.50 (LR/w, is the length of time-mean recirculation zone behind the body). The difference in length of the top and bottom recirculation zone shows the improper averaging space (spanwise). The streamlines centers were identified as local maxima of vorticity and saddle points were identified as minimal of vorticity (in all Figs. 3a to 3d, + represents the centre and × represents saddle) [9]. As studied by Lyn et al, [9], the two critical points may develop in the interior of the flow field (zero velocity or stagnation). The centres characterized by closed streamlines and saddles characterized by the intersection of streamlines. In the present case, near the cylinder, center of the large vortex lie at x'/w =1.12, y'/w = 0.35; and the saddle point was located at (x'/w, y'/w = 1.12, 0.95; the location is measured from the centre of the cylinder as shown in Fig. 3(a). Furthermore, the centers of the upper and lower bubbles located at approximately (x'/w =0.12, y'/w =0.6 and x'/w =0.08, y'/w = -0.62) and saddle point is located at (lower-side, x'/w = 0.6 and y'/w = -0.52), respectively.
The contours of normalized streamwise velocity (Ux/Uo) are shown in Fig. 3(c) which shows the negative values of velocity in the flow separation regions. Flow separates at the sharp leading edges of the square cylinder and the wake region is generated in the downstream region. The maximum value of velocity on the wake centreline is found to be Ux/Uo= 0.875. Zero velocity is clearly observed at the center and saddles point (as marked in Fig. 3(c). The contours of k normalized by Uo 2 for Case-I are shown in Fig.  3(d) and the highest value of k/Uo 2 =0.16, which occurs in the separating shear layers of the cylinder (top side, observed at x'/w =1.1 and y'/w =1.6). A saddle point is present near to this peak value. Moreover at the bottom edge of cylinder second-peak of k/Uo 2 =0.15 is observed at x'/w =0.08, y'/w = -0.62. A center and a saddle point are present near to this second-peak value. The value of k is reduced in the wake of the cylinder and at the streamline center of the large vortex (at x'/w =1.12, y'/w =0.35) k/Uo 2 =0.03. It shows that SST kω under-predict the value of Turbulent Kinetic Energy (TKE) in the wake region.
The normalized velocity (Ux/Uo) profiles for a cylinder in the wall-normal direction are shown in Fig. 4(a). The solid line shows the results obtained from the present SST k-ω model. Some deviations are observed in the near-wall region which shows the acceleration of velocity near the wall. In the far-field, the SST k-ω model shows good results when comparing it with previous results of Lyn et al. [9], Sohankar [18] and Nazari [29].
The distribution of streamwise normalized velocity (Ux/Uo) in the wake region of the cylinder (at y/w = 10.0) is presented in Fig. 4(b). It can be observed that Nazari et al. [29] over predicted these velocities in the far wake region. Moreover, the RSM results from  [9] data in the far wake region. Similarly, in the stream direction, the SST k-ω model does not show a descent match in the wake region of the cylinder and is unable to predict the data correctly. This may be due to an inappropriate size of the grid used in the wake region (see Fig. 1(b), or cause of SST k-ω URANS model that overpredicts the velocity in the wake of a single square cylinder. Whereas the present data shows approximate fit with the Nazari et al. [27] and Sohankar et al., [18] result in near wall region. The maximum deviation is observed between 11.8 ≤ x/w ≤ 13.8, as shown in Fig. 4(b). However, the present model show deviation when compared with the experiments of Lyn et al. [9]. Thus it is observed that the SST k-ω model fails to predict the accurate velocity fluctuations in the far wake region of a square cylinder.
The distribution of k/Uo 2 for Case-I is shown in Fig. 5. The SST k-ω predicts very low values of turbulent kinetic energy. Rodi [15] has observed similar distributions for k which show that RANS in general under-predicts k values. It is known that the k value is higher in separated flows which are not captured by RANS models accurately. Whereas, LES [15] is a capable technique to predict flow over bluff bodies as it shows good agreement with the experimental results of Lyn et al. [9] (Fig. 5). It is noted that the numerical predictions obtained by URANS simulations using SST k-ω model for flow over a square cylinder is not in agreement with experimental results. This is because such flow cases are inherently unsteady at high Reynolds number and the SST k-ω model is unable to predict these 3-D instabilities accurately.

Three Square Cylinders in Tandem (Case-II)
This section presents the results of the flow over three square cylinders placed in tandem named as Case-II. Firstly, the flow field from URANS is investigated and then the flow field from LES is analyzed. profiles along the centreline (at y'/w =0), Case-I Table 4(a) gives a summary of results obtained from the URANS simulations. The Case-II is compared with the case of flow over the two square cylinders investigated by Nazari et al. [29]. As the flow field is not exactly the same and there is a difference in Reynold's number (present, Re=22000, and Nazari et al. Re=16000). Thus this comparison is just used here to give a basic idea about the values of St and CD for a different number of square cylinders arranged in tandem. Whereas the present case is a novel exercise to simulate the flow around three square cylinders arranged in tandem. The value of St = 0.140 for both first and second cylinders, as reported by Nazari et al. [29]. The St related to vortex shedding in the two cases is different. However, from the present URANS simulations using the SST k-ω URANS model, the value of St from fine-grid are estimated to be around 0.1306 for the first cylinder (C-1) and 0.1299 for second cylinder (C-2), whereas it is 0.1004 for third cylinder (C-3, Table 4(b). The St is reduced for C-3 as it lies in the wake of the first two cylinders which causes the reduction in the vortex shedding from the last cylinder (C-3).

Flow Field from URANS
For the tandem arrangement of square cylinders Fig  6(a), the vorticity contour in the wake of the downstream cylinders shows qualitative differences with respect to that in the wake of the upstream cylinder mainly due to the difference in the nature of the flow that impinges upon them. The vortex peak coincides within the relevant streamline centers ( Fig.  6(b) in the wake of the square cylinder C-1. Whereas, streamlines saddles (×) in the wake correspond to zero vorticity (blue-contours, Fig. 6(a)). The velocity magnitude is also zero in all vortex cores (streamline centers) and at streamlines saddles points (see Figs. 7(d, e, f). It is mentioned by [18] that a region with lower pressure values is also defined as a vortex core (or center +, Fig. 6(c)). Thus the low-pressure zones in Fig. 6(c) are assumed to be the core of the vortex present in that region.
Moreover, Fig. 6(d) depicts the spatial evolution of the Turbulent Kinetic Energy (TKE) at a selected time instant. The contour shows that the generation of TKE is related to the shear in the flow arising from the cylinders. It is clear that the generation of TKE is greater in the regions where the velocity gradients are high, particularly in the near-wake (second-peak of TKE) of C-1 and the top side where the shear layer is separating (red-zone, first-peak of TKE). The values of the TKE around the cylinder C-2 and C-3 become small compared with those around the square cylinder C-1 (Fig. 6(d)). Near the C-1, the first saddle point is observed to be very close to the TKE second-peak, (see Fig. 6(d)) and it is further elaborated in the detail in Fig. 7(g).
The instantaneous pathlines for Case-II are shown in Fig. 7 (a. C-1, b. C-2, & c. C-3). Similar vortical structures are observed on C-1 in the tandem arrangement when comparing it with Case-I (singlecylinder of Lyn et al. [9]). Further analysis shows that as the flow is proceeding towards the last cylinder the Numerical; Nazari et al. [29] First cylinder (C-1) 16000 1.9299 Second cylinder (C-2) 0.8640 vortical structures are kept vanishing. Streamlines centers and saddle points are clearly observed for the streamlines around each cylinder. In the wake of C-1, a large rotating vortex appeares that strikes the second cylinder and the new vortex is generated in the wake of C-2 ( Fig. 7 (a, b, & c)). Thus the flow over C-2 is dependent on these vortical structures. Moreover, the aerodynamic forces are also depending on the vortex shedding. It can be seen from Table 4(a) that the values of drag coefficients (CD) for C-2 and C-3 are 0.583 and 0.249 and these are much smaller than that of C-1 (CD=2.372). So overall a decrease in the drag forces is observed for the downstream cylinders that lie in the wake of upstream cylinders.
The distribution of streamwise velocity Ux/Uo obtained from URANS simulation using SST k-ω model is presented in Fig. 7 (d, e, & f). This instantaneous velocity distribution shows quite asymmetrical behavior. Higher values of velocities are observed at the top and bottom sides of each cylinder. Moreover, negative values are reported in flow separation regions and the wake regions of cylinders.
Whereas zero velocities are clearly captured near centre and saddles of streamlines. The overall gradual decrement in velocity is monitored as flow moves from the cylinder C-1 to C-3. The turbulence is further demonstrated by the distribution of normalized k contours in Fig. 7(g, h, AND i). The highest value of Fig. 7: Streamlines, Streamwise Velocity (UX/UO), and Turbulent Kinetic Energy (k/UO 2 ) Contours, Case-II-URANS, (at z/w=0.5) (a, d, g) C-1, (b, e, h) C-2, (c, f, i) C-3; (Here, + represent centre, and × epresent saddle) k/Uo 2 =0.15 is reported at the top side of C-1, as the flow is accelerating in this region and shear stresses are higher in this zone. A streamline centre is also present near this peak-value whereas a saddle point is almost coinciding with this first-peak of TKE ( Fig.  7(g)). Further, a second-peak value of k/Uo 2 =0.11 is detected at x'/w =1.4, y'/w = -0.85 and a visible saddle point is also present near to this second-peak (Fig. 6(d) and Fig. 7(g)). This observation is not in accordance with the [9] study. As Lyn et al. observed opposite behaviour of streamline centre and saddle points (for details see the [9] study and next section 5.2.2 LES). Besides, low values are observed in the wake of C-1. Moreover, as the flow passes the first cylinder and reaches the second cylinder, k decreases and is weak around C-2. Furthermore, the vortex shedding and turbulence becomes weaker as the flow approaches the cylinder C-3 and a low value of about 0.03 is observed in the wake of C-3. For all cylinders, the k is higher at the leading edges (top and bottom corners) of cylinders C-1, C-2, and C-3, whereas the value of k decreased in wake regions of cylinders.

Flow Field from LES
Instantaneous flow structures obtained from LES are studied using the second invariant of velocity gradient tensor (Q). The iso-surfaces of Q distribution {Q = 0.5 (u #,# u $,$ − u #$ u $,# )} are presented in Figs. 8 (a, b, c and d) at four instants. The flow shows the unsteady character and the wake oscillates in time and space. The complex flow field causes the variation in drag and lifts experienced by the cylinders. Large Q values (Q > 0) indicate the strong vortical flows. Periodic vortex shedding is clearly observed at the four instants of time. Some grid resolution errors can be observed around the cylinder due to highly skewed grid cells present in the near wake and wall-normal directions. Thus the size of the grid can be modified to improve the ability of solver for capturing good quality of the flow structures (the study of using different grid sizes for LES will be performed in future work). However, presently the DSM-LES method is found to be efficient in capturing the different sizes of the vortical structures present in the surroundings of the three square cylinders.
The plots of pathlines for Case-II (LES) are shown in Fig. 9 (a, b and c). Vortices shed at the upper and lower faces of the cylinders and the flow separates from the leading edges of the cylinders. A total of four vortices can be seen in the vicinity of each cylinder. These include two vortices one at the top, second at bottom of the cylinder and a pair of counter-rotating vortices in the near wake of the cylinders. The size of vortices at the top and bottom of the cylinder reduces as the flow sweeps over the row of cylinders. The size of vortices in the wake region for the middle cylinder is larger than the first and the third cylinder. In the wake region, the length of the recirculation zone (LR/w) for C-1 is 0.78, for C-2 it is 0.98 and for C-3, it is 0.80. From the analysis of pathlines ( Fig. 9 (a, b, and c)), it is observed that center and saddle points are present in wake regions, and at the top and bottom sides of each cylinder. The value of the velocity at these centre and saddles points is observed to be zero from the velocity contours in Fig. 9 (d, e, and f). The negative velocities are detected in the wake regions of cylinders placed in the tandem arrangement (as shown in Fig. 9 Fig. 10 (a, b, and c). As the flow moves from one cylinder to the next cylinder, the decrement in values of u v are perceived. The maximum values of u v for C-1 is -0.16 at (x'/w =1.8, y'/w =0.5), for C-2 it is 0.11 at (x'/w =0.5, y'/w = -0.3), for C-3 is 0.06 at (x'/w =0.2, y'/w = -0.6) in the wake region. Similar trends are observed by Sohankar [18] in the investigation of flow over a single square case at Reynolds number of 2. Observed to be away from the centre point (vorticity peak) for all cylinders (in accordance with the experimental study of Lyn et al. [9]). The locations of high shear stresses and high normal stresses coincide with underlying events of the vortex shedding, and the significant turbulence production originating from the separating shear layers can be observed (see Fig. 11(ato-f)).
The contours of stream-wise turbulent normal stresses u u for the flow passing the square cylinders placed in tandem are presented in Fig. 10 (d, e and f). Fairly symmetrical structures are observed around the centreline of each cylinder. High levels of u u are found at the sides of cylinders. Furthermore, the values of normal stresses u u are much higher than that of the shear stresses u v . Moreover, as mentioned earlier that turbulent kinetic energy (k) is defined as k = (u u + v v + w w ). The distribution of k is shown in Figs. 10 (g, h and i) for square cylinders in tandem. It is observed that k is highly dependent on u u components. Normal stress (F F ) surrounding the last cylinder is small in magnitude as compared to the first and the second cylinder. Similarly, for C-1 the value of k is around 0.36 located at x/w=12.20 because the wake vortices from the upper and lower shear layers pass through almost the same position close to the centreline. Here the distribution of u u (Fig.  10(d)) shows that the upper and lower vortices intersect at the centre line.
The maximum value of TKE is found in the wake of cylinder C-1 ( Fig. 10(g)). It means that evident vortex shedding is initiated from the upstream cylinder C-1. For C-2 TKE is found to be 0.21 at x/w=18.60 and for C-3 the value is 0.12 at x/w=25.80. Thus k is gradually reduced from C-1 to C-3. Moreover the peak values of k and u u appeared near to the centres of the streamlines vortices, thus it contradicts the present URANS results of TKE whereas it is according to Lyn et al. [9] work. On the other hand, the centre and saddle points surrounding the cylinder C-2 (top and bottom, Fig. 10(h)) are again showing that these points are near to the peak-TKE. Also the last cylinder C-3 peak-TKE value of 0.12 is overlapping a saddle point ( Fig. 10(i)).
An in-depth analysis of turbulent production due to shear stresses and normal stresses is detailed here. The production of turbulent kinetic energy due to shear stress Pks and normal stress Pkn can be evaluated Pkn are shown in Fig. 11 (a, c, e) and Fig. 11 Fig. 11 (b). Similarly for C-2 peak of Pks lies near the top and bottom of the cylinder and smaller values are present in the wake of the cylinder (Fig. 11  c). As compared to C-1 and C-2 the Pks values are significantly small around C-3 ( Fig. 11(e)).
In Figs. 11(b, d, f), the negative Pkn values (dotted line) are more prominent; this arises in the regions where v M u [9]. The negative peak value of the production of turbulent kinetic energy due to normal stress (Pkn= -0.15) exactly coincides with the saddle point in the wake of cylinder C-1 ( Fig. 11(b)). Whereas negative peak values are also detected near cylinder C-1 top (Pkn= -0.4) and bottom (Pkn= -0.45) regions and vorticity-centers (+) are also present near to it ( Fig. 11(b)), thus it is in agreement with the study of Lyn et al. [9]. Around all cylinders, the regions of negative Pkn are also consistently present in vortexpeak zones [9]. Lastly, it is clearly observed that the Pks and Pkn contours are elongated for C-1 as compared with the last two cylinders. Thus it confirms the vortex shedding is initiated from first cylinder C-1 and these vortexes further impinge on the second cylinder C-2 and lose its turbulence before striking the cylinder C-3. That is why the turbulence around C-3 decays, and very small values of Pks and Pkn are observed around C-3 (Figs. 11(e and f)). An overall analysis of turbulence production shows that the main features of these complex flows can be predicted reasonably well by the present-LES Dynamic Smagorinsky model (DSM).

Comparison of URANS and LES for Case-II
In this section, the case of flow over three square cylinders placed in the tandem arrangement is compared with the experimental and numerical data of Case-I. The present Case-II is a novel investigation and this has been never done before. Thus, these comparison graphs (shown below in Fig. 12, 13 and 14) will only approximate the results related to the flow over the first cylinder (C-1) with the case studies of single square cylinder; whereas, the graphs shown for cylinders C-2 and C-3 are only comparison of present URANS and LES models.
The streamwise velocity distribution at three different locations for flow over the square cylinders placed in tandem is presented in the form of a graph in Fig. 12.
The LES predictions and URANS simulation using the ST k-ω model are compared with the experimental data of Lyn et al. [9]. It is observed that the velocity distributions at the first cylinder (C-1) in the tandem arrangement are very similar in trend and magnitude of the flow over a single square cylinder. The minimum velocities are located at y'/w=0.5 from the center of the first cylinder with values -0.24, -0.20 and -0.15 for Lyn et al. [9], LES, and URANS, respectively. It is also found that separated shear layers do not reattach at any point on the upper and lower edges of cylinder C-1 as observed in Fig. 9 (a) ( Fig.  12(a), is in agreement with Rodi [15]). Here LES calculations are in agreement with the experiment of [9]; whereas SST k-ω model shows deviations.
Overall LES technique examined the appropriate values when comparing it with experimental results of Lyn et al. [9] with a percentage deviation of approximate 17% in the near-wall region and with a 10% deviation in the far wake region. Whereas URANS showed deviations of 37.5% in the near-wall region and 10% in far wake regions. Moreover, from the plots of cylinders C-2 and C-3, the URANS observations are in accordance with the LES (both shows an increase in velocity in the streamwise directions). Streamwise velocity is plotted along the centreline for Case II in Fig. 13. From these graphs, it is accessed that in comparison to Lyn et al. [9], LES is accurately predicting the velocity fluctuations in the near-wall region and suitable readings are also captured in the far wake of the first cylinder ( Fig. 13(a)). Besides that URANS overpredicted the values of Ux/Uo. An inappropriate match is detected in LES and URANS techniques from the graph of C-2. However, in graph C-3 LES intersect URANS data and nearly predict similar values for velocities. The graphs for normalized turbulent kinetic energy (k/Uo 2 ) fluctuations over the three cylinders are shown in Fig.  14. Maximum values of k/Uo 2 are observed in the wake of the first cylinder whereas minimum values of k/Uo 2 are monitored in the wake of the third cylinder, and this is the evidence that low turbulence is present around the last cylinder (C-3). In comparison to experiment of Lyn et al. [9] and LES results of Rodi [15] the graph of k/Uo 2 for the C-1 shows that SST kω model is unable to capture these fluctuations compared to the previous results as mentioned in Figs.12a and 13a, whereas LES, highly over predict these values. As Rodi [15] has investigated flow over a single cylinder, and our Case-II the C-1 is aligned with two more cylinders close to each other (tandem arrangement), thus the TKE is not similar to single square cylinder case. Moreover, the differences in LES and URANS by this graphical presentation of k/Uo 2 are much higher for the second and the third cylinder. That confirms the differences present in the simulation results captured by these two techniques. These results also show that the complex flow-field is very difficult to be realistically simulated by the SST k-ω URANS turbulence model.

CONCLUSIONS
The flow field over the square cylinders is very complex due to inherent flow separation unsteadiness and anisotropy. There are distinct coherent structures present in the separated flow regions which are undulating in time and space. The flow field over a single square cylinder and three square cylinders placed in a tandem arrangement are investigated using different numerical models. LES and URANS based models are tested and compared with the experimental work. The Reynolds number used in the present work is 22,000.
The URANS based turbulence model is found to predict the frequencies of vortex shedding that are not in order with the Lyn et al. [9]. As the numbers of cylinders are increased the complexity of the flow field increases and the SST k-ω turbulence model, due to its limitations, has failed to capture the true physics of the problem. The model is also failed in predicting the true levels of turbulent kinetic energy (k). The SST k-ω URANS-model results are not entirely satisfactory and do not provide complete turbulence analysis and periodic motion is strongly unpredicted. One reason for under predicting the values of TKE could be the use of z/w=1.0 in the spanwise direction and this issue will be considered in future-work Rodi [15].
However, the Dynamic Smagorinsky LES-model (DSM) shows improved results. LES based model predicts the 3-D fluctuations in a proper manner and in general gives better details of the flow. Present LES confirms that the production of turbulence by normal stresses is higher than the production of turbulent