# Numerical Study on Flow around Four Square-Arranged Cylinders at Low Reynolds Numbers

Authors: Yang-yang Gao; Chang-shan Yin; Hao-qiang Zhang; Kang Yang; Xi-zeng Zhao; Zhilin Sun
Journal: Mathematical Problems in Engineering
 (2017)

---

## Abstract

In this paper, numerical simulations of flow past four square-arranged cylinders are carried out at different spacing ratios (1.5≤L/D≤5.0; L is the center to center distance; D is the cylinder diameter) and Reynolds numbers (100≤Re≤1000). The effects of spacing ratio and Reynolds number on the wake flow characteristics are investigated, such as the instantaneous vorticity contours, force coefficients, and vortex shedding frequencies. The results show that the flow characteristics behind the four-cylinder cases are significantly affected by the spacing ratios and Reynolds numbers. At the same spacing ratio, the transformation of flow pattern is advanced quickly with increasing of Reynolds numbers, the values of force coefficients are correspondingly fluctuated with large amplitude, and the vortex shedding frequency is increased significantly with Re.


---

## Body

## 1. Introduction

Investigation of unsteady viscous flow past a single or array of circular cylinders is of great interest for many ocean and offshore engineering applications. Compared to an isolated single cylinder, it is more difficult to predict flow around multicylinder configurations, depending on the geometric parameters such as spacing ratio between the cylindersL/D (L is the center-center distance between cylinders, and D is the cylinder diameter) and the incidence angle α relative to the free-stream flow. The flow behaviors behind multicylinders are much more complex due to the complicated dynamic interaction between the shear layers, shed vortices, and Karman vortex streets in the wake flow. Most of the researches are focused on flow around an isolated circular cylinder or two cylinders, as pointed in the reviews by Zdravkovich [1], Sumner [2], and Bearman [3]. However, few investigations are performed to reveal the flow dynamics of a four-square-arranged cylinder configuration, which is essential for various marine structures as a fundamental arrangement. Four square-arranged cylinders could be found in many offshore engineering applications, such the pile foundation system of offshore platforms and offshore wind turbine, risers, and pipelines. When the four cylinders are in close proximity configuration, flow-induced vibration may arise, which would increase the risk of fatigue failure in engineering. Therefore, the numerical simulation was carried out to investigate the flow characteristics of flow past four square-arranged cylinders.The earliest experimental investigations on flow past four square-arranged cylinders were reported by Sayers [4, 5]. The force coefficients and vortex shedding frequency for each cylinder were investigated at 1.1≤L/D≤5.0 and Re=3.0×104 (Reynolds number is defined as Re=UD/v, D is the diameter of cylinder, U is the free-stream velocity, and v is the kinematic viscosity) in an open-jet wind tunnel. The nondimensional vortex shedding frequency St (defined St=fD/U; f is the vortex shedding frequency) was found to be varied largely across wake at L/D<4.0 and equal to that for a single isolated cylinder when L/D≥4.0. Soon after that, Lam and Lo [6] categorized the flow characteristics behind four square-arranged cylinders into three flow patterns varied with L/D ranging from 1.28 to 5.96 at Re=2100: (i) the shear layers of upstream cylinder shielded the downstream cylinder; (ii) the generated free shear layers from upstream cylinder reattached onto the downstream cylinder; (iii) the vortices shed from upstream cylinder and impinged the downstream cylinder. Lam and Fang [7] experimentally studied force coefficients and the effects of flow interference between four cylinders on the mean pressure distribution while L/D ranges from 1.26 to 5.80 with H/D=28.4 and Re=12.8×103. They confirmed that L/D=2.7 is a critical spacing ratio and the forces characteristics vary largely with incidence angles (0°  ≤ α  ≤ 45°) and L/D. Similarly, Lam et al. [8] experimentally measured the mean and fluctuating forces and Strouhal numbers on each cylinder at spacing ratio L/D in the range 1.69 to 3.83, different angles of incidence α ranging from 0° to 180° at a 15° interval, and Re=41.0×103. As a result, three types of flow pattern were observed depending on L/D: (i) shear layers around cylinder oscillate freely and vortex forms at large L/D, (ii) shear layers reattach on or shield the downstream cylinder, and (iii) narrow gap flows form at small L/D. Lam and Zou [9, 10] summarized several distinct flow patterns at subcritical Reynolds numbers (11000≤Re≤20000) and concluded that the turbulence flow patterns are affected by the spacing ratio and Re by using LDA (laser Doppler anemometry) and PIV (particle imaging velocimetry) to measure the turbulent flows around four cylinders in an in-line square configuration with different spacing ratios of 1.5, 2.5, 3.5, and 5.0. Lam et al. [11] carried an experimental investigation on flow past four cylinders in a square configuration with a spacing ratio of 4.0 and at Re=200 by using PIV and LIF (laser induced fluorescence) flow visualization technique. Concerning α, they observed three basic flow patterns including the formation of a jet flow which could lead to strong flow-induced vibration. Recently, Wang et al. [12] further investigated the vortex shedding characteristics using PIV at various L/D from 2.0 to 5.0, Re=8000, and inclination angle ranging from 0° to 45°. The fluctuating drag and lift coefficients on four cylinders were measured and the relation between the flow patterns and forces on the cylinders was further studied. By using LIF and PIV technology, Zou and Lin [13] captured the biased flow pattern and classical bistable phenomenon for flow past four square-arranged cylinders at Re=200, L/D=1.6, and H/D=16 (H is the length of cylinder immersed in water), which is seldom investigated by researchers.Apart from the above experimental studies, numerical simulations on four-cylinder array have also been investigated at small Reynolds numbers, except for some LES study [9, 10, 14] at Re=1.5×104 and the FEM study by Zhao and Cheng [15] at Re=103–2×104. Farrant et al. [16] employed the cell boundary element method (CBEM) to investigate the two-dimensional incompressible flow around four-cylinder array at Re=200 and inclination angle α = 0° and 45°. The in-phase and anti-phase vortex shedding behavior were observed in their study. Lam et al. [17] carried out a two-dimensional simulation on cross-flow past four square-arranged cylinders at Re=100 and 200. The results of their study showed that, depending on the spacing ratio L/D, three distinct flow patterns were observed, and the alternation of forces and pressure distribution on cylinders when different flow patterns transform were also discussed. Han et al. [18] used spectral element method (SEM) to simulate a two-dimensional laminar flow (Re=200) with incidence angle α = 0° and 45°. Three distinct flow patterns exist with variation of L/D in each incidence angle, and the force coefficients and Strouhal numbers changed significantly when one kind of flow patterns switched to another. Abbasi et al. [19] further performed 2D numerical investigation of flow past four square cylinders in an in-line square configuration using LBM (lattice Boltzmann method). Four distinct flow patterns were observed at different Reynolds numbers (60≤Re≤175) when spacing ratio was set at 2.0, 4.0, and 7.0, respectively. They also noted that the variation of Re was contributed much to the transformation of flow pattern and the characteristics of wake flow when spacing ratio is small (L/D=2.0). Tong et al. [20] carried out a three-dimensional simulation on flow past four square-arranged cylinders at L/D=2.0 and 100≤Re≤500. Four wake flow regimes were observed at different Re, and they found that the root-mean-square (RMS) values of dynamical parameters, length of wake flow, and phase angle of lift coefficient on downstream cylinders were all changed significantly when the flow transited from one regime to another. Generally, the published works discussed above are summarized in Table 1.Table 1
The experimental results and calculated results for the flow around four cylinders.
Investigators
Re
L
/
D
α
Technique
Experiments
Sayers [4, 5]
3.0 × 104
1.1–5.0
0°–180°
Pressure method
Lam and Lo [6]
2.1 × 103
1.28–5.96
0°–45°
LIF, PIV
Lam and Fang [7]
1.28 × 104
1.26–5.80
0°–45°
Pressure method
Lam et al. [11]
200
4.0
0°–45°
LIF, PIV
Lam et al. [8]
4.1 × 104
1.69–3.83
0°–180°
PELC
Lam and Zou [9]
(1.128–1.982) × 104
1.5–5.0
LDA, PIV
Lam and Zou [10]
(1.1 –2.0) × 104
1.5–5.0
LDA, PIV
Zou and Lin [13]
200
1.6
0°
LIF, PIV
Wang et al. [12]
8000
2.0–5.0
0°–45°
PIV, PELC
2D simulations
Zhao and Cheng [15]
1.0 × 103–2.0 × 104
3.0
0°–45°
FEM
Farrant et al. [16]
200
3.0, 5.0
0°, 45°
CBEM
Lam et al. [17]
100, 200
1.6–5.0
FVM
Han et al. [18]
200
1.5–4.0
0°, 45°
SEM
Abbasi et al. [19]
60–175
2.0, 4.0, 70
LBM
3D simulations
Lam and Zou [9]
1.5 × 104
1.5, 3.5
LES
Lam et al. [17]
200
4.0
FVM
Zou et al. [14]
1.5 × 104
1.5, 3.5
LES
Lam and Zou [10]
1.5 × 104
1.5, 3.5
LES
Tong et al. [20]
100–500
2.0
FVM
LIF: laser induced fluorescence; PIV: particle image velocimetry; LDA: laser Doppler anemometry; PELC: piezoelectric load cell; FEM: finite element method; CBEM: cell boundary element method; FVM: finite volume method; SEM: spectral element method; LBM: lattice Boltzmann method; LES: large eddy simulation.It is important to note that most of the mentioned previous experimental and numerical studies on flow past a four-cylinder configuration are carried out at relatively high subcritical Reynolds numbers or variation with spacing ratio at a constant Reynolds number. However, few studies about the combined effects of Reynolds number and spacing ratioL/D on wake flow characteristics are presented. In this paper, two-dimensional simulations on flow past four circular cylinders in an in-line square configuration are carried out using the finite volume method (FVM). The spacing ratios L/D vary from 1.5 to 5.0 at a 0.5 interval and Reynolds numbers are varied from 100 to 1000. Concerning the laminar flow (Re≤200) and subcritical flow (Re≥300) in this study, different proper calculation models are adopted for them, respectively. In the present work, ANSYS Fluent software was used to conduct simulation work. The objective of this work is to provide further insight into the combined effects of spacing ratio L/D and Reynolds number on the wake flow structures and dynamic responses of four square-arranged cylinders. The complicated interactions between the wake flow and the multiple cylinders are also discussed.


## 2. Numerical Method

### 2.1. Governing Equations and Models

In this numerical simulation, laminar model was employed forRe=100 and 200, and the nondimensional tensorial forms of the Navier-Stokes governing equations for the incompressible viscous fluid flow can be expressed as(1)∂ui∂t+uj∂ui∂xj-1Re∂2ui∂xj∂xj+∂p∂xi=0,∂ui∂xi=0.ui  (i=1,2) are the velocities in the 2-dimensional Cartesian coordinates (xi or xj (i,j=1,2)), p is the pressure, t is the time, and Re is the Reynolds number. For the cases of Re≤200, simulation could be conducted by directly solving (1).WhenRe=300, 500, and 1000, the standard k-ω model derived from Wilcox [21] is employed and the nondimensional tensorial forms of the time-averaged Navier-Stokes governing equations for the viscous fluid flow can be expressed as (2)∂ρ∂t+∂∂xjρuj=0,(3)∂∂tρui+∂∂xjρujui=-∂p∂xi+∂∂xjμ∂ui∂xj-ρui′uj′¯.-
ρ
u
i
′
u
j
′
¯ in (3) is Reynolds stress and it could be simplified as(4)-ρui′uj′¯=μt∂ui¯∂xj+∂uj¯∂xi-23δijρkwhen the Boussinesq hypothesis [22] is considered. However, more variables are introduced in (3), so, in order to make (2) and (3) closed for solution, two equations are introduced, namely, the transport equation of k and ω: (5)∂ρk∂t+∂∂xjρujk-μ+σ∗μt∂k∂xj=τij∂ui∂xj-β∗ρωk,∂ρω∂t+∂∂xjρujω-μ+σμt∂ω∂xj=αωkτij∂ui∂xj-βρω2.The definitions of related parameters in (3) and (5) are given as (6)τij=2μtSij-13∂uk∂xkδij-23ρkδij,Sij=12∂ui∂xj+∂uj∂ui,μt=α∗ρkωwith the corresponding explanation given as follows: Sij, mean strain rate tensor; τij, Reynolds stress tensor. And the definitions of the other parameters are presented in the following:ρ, fluid density; μ, the molecular viscosity; δij, the Kronecker delta; μt, the eddy viscosity; k, the turbulent kinetic energy; ω, the specific dissipation rate. The values of some factors in (5) and formula (6) can be referenced to the turbulence model proposed by Wilcox [21].When-ρui′uj′¯ is solved using (5), and (2) and (3) could be solved when the solution of -ρui′uj′¯ is taken into (2) and (3). The N-S equations combined with turbulent model are finally solved.In addition, Stringer et al. [23] suggested the turbulence intensity I to be 2.5% for the numerical simulations of flow past a single circular cylinder at different Reynolds numbers. Hereby, some adopted parameters are as follows:(7)k=32U∞I2,ω=k3/2kCμl,l=0.07D,kwall=0,ωwall=60μρβy12.


### 2.2. Computational Domain and Boundary Conditions

Figure1 shows the schematic of the geometry and mesh of flow past four square-arranged circular cylinders. The computational domain is 60D×40D (D is the diameter of cylinder) and the coordinate origin (0,0) is set at the center position of the four-cylinder arrangement. The inlet of the computational domain is placed 20D upstream from the coordinate origin and the outlet is placed 40D downstream from the coordinate origin. The upper and lower boundary are placed 20D away from the coordinate origin, with symmetry condition, which is also adopted in the numerical simulation of flow past four square-arranged circular cylinders by Han et al. [18], Tong et al. [20], and Ji et al. [24]. The physical meaning of the symmetry boundary condition has been presented in Table 2, the normal gradient of stream velocity is zero, and lateral velocity is zero. The four-cylinder array could not be affected by the side wall, which could ensure that the cylinders are normal to uniform flow. The dimensions of domain in this study can be referred to Stringer et al. [23]. In addition, no-slip conditions are adopted for the four cylinders’ surface. The details for the boundary conditions are shown in Table 2. The velocities of uniform flow at different Reynolds numbers and the corresponding calculation models are summarized in Table 3.Table 2
Boundary conditions for computational domain.
Inlet
U
=
U
∞
,  V=0
Outflow
∂
U
/
∂
x
=
0,  ∂V/∂x=0
Symmetry
V
=
0,  ∂U/∂y=0
No-slip
U
=
0,  V=0Table 3
Calculation models at different Reynolds numbers.
Re
U
∞ (m/s)
Calculation model
100
0.005
Laminar
200
0.010
Laminar
300
0.015
Standardk-w
500
0.025
Standardk-w
1000
0.050
Standardk-wFigure 1
The schematic of the computational domain (a), and the typical local mesh near four cylinders (b). The top-right corner of the second image is the zoomed in view close to the cylinder surface.
(a)
(b)In this paper, laminar flow was simulated for the case ofRe=100. At Re=200, numerous researchers carried out 2D laminar flow simulations for multicylinder arrays, such as Farrant et al. [16], Lam et al. [17], Harichandan and Roy [25], and Han et al. [18]. The values of C¯D and St computed from laminar flow are close to the results computed from turbulent flow model [26]. The difference for C¯D between Vu et al. [26] and current work is below 0.6%, and the difference for St is below 3%. Therefore, the effect of turbulence could be ignored in the present study. However, to get accuracy results, turbulent calculation model was used at Re≥300.In this study, finite volume method (FVM) was employed to discretize these equations. The well-known SIMPLEC algorithm was used to deal with the coupling between the pressure and the velocity fields, which could be helpful to speed up the convergence of equation and ensure agreeable accuracy during the calculation [27]. As for the discretization of the convective or pressure terms in the conservations, it was accomplished through a second-order upwind differencing scheme, which has been adopted by Lam et al. [17]. The second-order upwind differencing scheme was used in this study due to its high stability and veracity compared to the first-order upwind scheme (see Shyy et al. [28]; Ni et al. [29]). And it is the same for transport equation of k and ω in terms of the discretization scheme. The second-order implicit forward discretization was adopted for the time derivative term in order to reduce the convergence time. The continuity equation is more difficult to converge than other equations; therefore the criteria for convergence were set in a manner that residual is 10−4 for continuity equation and 10−5 for other equations, respectively. The schematic of calculation procedure based on the SIMPLEC algorithm is presented in Figure 2.Figure 2
The schematic of calculation procedure.


## 2.1. Governing Equations and Models

In this numerical simulation, laminar model was employed forRe=100 and 200, and the nondimensional tensorial forms of the Navier-Stokes governing equations for the incompressible viscous fluid flow can be expressed as(1)∂ui∂t+uj∂ui∂xj-1Re∂2ui∂xj∂xj+∂p∂xi=0,∂ui∂xi=0.ui  (i=1,2) are the velocities in the 2-dimensional Cartesian coordinates (xi or xj (i,j=1,2)), p is the pressure, t is the time, and Re is the Reynolds number. For the cases of Re≤200, simulation could be conducted by directly solving (1).WhenRe=300, 500, and 1000, the standard k-ω model derived from Wilcox [21] is employed and the nondimensional tensorial forms of the time-averaged Navier-Stokes governing equations for the viscous fluid flow can be expressed as (2)∂ρ∂t+∂∂xjρuj=0,(3)∂∂tρui+∂∂xjρujui=-∂p∂xi+∂∂xjμ∂ui∂xj-ρui′uj′¯.-
ρ
u
i
′
u
j
′
¯ in (3) is Reynolds stress and it could be simplified as(4)-ρui′uj′¯=μt∂ui¯∂xj+∂uj¯∂xi-23δijρkwhen the Boussinesq hypothesis [22] is considered. However, more variables are introduced in (3), so, in order to make (2) and (3) closed for solution, two equations are introduced, namely, the transport equation of k and ω: (5)∂ρk∂t+∂∂xjρujk-μ+σ∗μt∂k∂xj=τij∂ui∂xj-β∗ρωk,∂ρω∂t+∂∂xjρujω-μ+σμt∂ω∂xj=αωkτij∂ui∂xj-βρω2.The definitions of related parameters in (3) and (5) are given as (6)τij=2μtSij-13∂uk∂xkδij-23ρkδij,Sij=12∂ui∂xj+∂uj∂ui,μt=α∗ρkωwith the corresponding explanation given as follows: Sij, mean strain rate tensor; τij, Reynolds stress tensor. And the definitions of the other parameters are presented in the following:ρ, fluid density; μ, the molecular viscosity; δij, the Kronecker delta; μt, the eddy viscosity; k, the turbulent kinetic energy; ω, the specific dissipation rate. The values of some factors in (5) and formula (6) can be referenced to the turbulence model proposed by Wilcox [21].When-ρui′uj′¯ is solved using (5), and (2) and (3) could be solved when the solution of -ρui′uj′¯ is taken into (2) and (3). The N-S equations combined with turbulent model are finally solved.In addition, Stringer et al. [23] suggested the turbulence intensity I to be 2.5% for the numerical simulations of flow past a single circular cylinder at different Reynolds numbers. Hereby, some adopted parameters are as follows:(7)k=32U∞I2,ω=k3/2kCμl,l=0.07D,kwall=0,ωwall=60μρβy12.


## 2.2. Computational Domain and Boundary Conditions

Figure1 shows the schematic of the geometry and mesh of flow past four square-arranged circular cylinders. The computational domain is 60D×40D (D is the diameter of cylinder) and the coordinate origin (0,0) is set at the center position of the four-cylinder arrangement. The inlet of the computational domain is placed 20D upstream from the coordinate origin and the outlet is placed 40D downstream from the coordinate origin. The upper and lower boundary are placed 20D away from the coordinate origin, with symmetry condition, which is also adopted in the numerical simulation of flow past four square-arranged circular cylinders by Han et al. [18], Tong et al. [20], and Ji et al. [24]. The physical meaning of the symmetry boundary condition has been presented in Table 2, the normal gradient of stream velocity is zero, and lateral velocity is zero. The four-cylinder array could not be affected by the side wall, which could ensure that the cylinders are normal to uniform flow. The dimensions of domain in this study can be referred to Stringer et al. [23]. In addition, no-slip conditions are adopted for the four cylinders’ surface. The details for the boundary conditions are shown in Table 2. The velocities of uniform flow at different Reynolds numbers and the corresponding calculation models are summarized in Table 3.Table 2
Boundary conditions for computational domain.
Inlet
U
=
U
∞
,  V=0
Outflow
∂
U
/
∂
x
=
0,  ∂V/∂x=0
Symmetry
V
=
0,  ∂U/∂y=0
No-slip
U
=
0,  V=0Table 3
Calculation models at different Reynolds numbers.
Re
U
∞ (m/s)
Calculation model
100
0.005
Laminar
200
0.010
Laminar
300
0.015
Standardk-w
500
0.025
Standardk-w
1000
0.050
Standardk-wFigure 1
The schematic of the computational domain (a), and the typical local mesh near four cylinders (b). The top-right corner of the second image is the zoomed in view close to the cylinder surface.
(a)
(b)In this paper, laminar flow was simulated for the case ofRe=100. At Re=200, numerous researchers carried out 2D laminar flow simulations for multicylinder arrays, such as Farrant et al. [16], Lam et al. [17], Harichandan and Roy [25], and Han et al. [18]. The values of C¯D and St computed from laminar flow are close to the results computed from turbulent flow model [26]. The difference for C¯D between Vu et al. [26] and current work is below 0.6%, and the difference for St is below 3%. Therefore, the effect of turbulence could be ignored in the present study. However, to get accuracy results, turbulent calculation model was used at Re≥300.In this study, finite volume method (FVM) was employed to discretize these equations. The well-known SIMPLEC algorithm was used to deal with the coupling between the pressure and the velocity fields, which could be helpful to speed up the convergence of equation and ensure agreeable accuracy during the calculation [27]. As for the discretization of the convective or pressure terms in the conservations, it was accomplished through a second-order upwind differencing scheme, which has been adopted by Lam et al. [17]. The second-order upwind differencing scheme was used in this study due to its high stability and veracity compared to the first-order upwind scheme (see Shyy et al. [28]; Ni et al. [29]). And it is the same for transport equation of k and ω in terms of the discretization scheme. The second-order implicit forward discretization was adopted for the time derivative term in order to reduce the convergence time. The continuity equation is more difficult to converge than other equations; therefore the criteria for convergence were set in a manner that residual is 10−4 for continuity equation and 10−5 for other equations, respectively. The schematic of calculation procedure based on the SIMPLEC algorithm is presented in Figure 2.Figure 2
The schematic of calculation procedure.


## 3. Results and Discussions

### 3.1. Validation Studies

In order to validate the accuracy of the present numerical model, two-dimensional simulation on flow past a single circular cylinder is carried out at different Reynolds numbersRe=200, 300, 500, and 1000. Figure 3 shows the schematic of the computational structured mesh and the local mesh near the cylinder surface. Four different structured finite element meshes based on the number of nodes on the cylinder surface are compared in Table 4. It can be seen that when the number of nodes on cylinder surface exceeds 140, the values of mean drag coefficient C¯D, root-mean-square drag, and lift coefficient CD′, Cl′ and Strouhal number St are almost stable. To reduce consumption of computational resources, the nodal point of 140 instead of 180 on the circumference of the cylinder surface was adopted in this study. The O-Grid technique is applied to refine the quality of mesh around the cylinder, where the minimum mesh size is 0.002D mm. The nondimensional time step (Δt=Ut/D) is 0.001 in all of the computational cases.Table 4
Grid independence test atRe=200.
Number of nodes on cylinder surface
Number of cells
Δ
t
C
¯
D
C
D
′
C
l
′
St
Grid 1
100
11328
0.001
1.330
0.030
0.466
0.193
Grid 2
140①
16048
0.001
1.337
0.031
0.474
0.191
Grid 3
180
20768
0.001
1.340
0.032
0.485
0.191
Grid 4
220
25488
0.001
1.340
0.032
0.486
0.192
①The final chosen calculation scheme.Figure 3
The view of structured mesh around the single circular cylinder. (a) Full view of the computation domain; (b) detail mesh close to the cylinder surface.
(a)
(b)Table5 shows the comparison of results for flow past a single cylinder from other numerical simulations at Re=200. For this Re, the present result of C¯D,CD′,Cl′, and St is consistent with the previous results obtained from Ding et al. [30], Harichandan and Roy [25], Lam et al. [17], and Mittal [31]. Compared with the result of Ding et al. [30], the value deviation is 0.82% for C¯D, 12.90% for CD′, 1.93% for Cl′, and 2.62% for St, respectively. When Reynolds number Re≥300, the comparison of results from different turbulent models is carried out. Previous numerical results for flow past single circular cylinder at different Re ranging from 300 to 1000 are presented in Table 6. Different simulation codes and grids may cause a bit disparity between other works [20, 23, 26, 31–35] for related parameters listed in Table 6. However, it can be seen from Table 6 that the results obtained from the Standardk-w model appear to be more consistent with the other numerical studies for the cases of Reynolds number Re=300, 500, and 1000.Table 5
Comparable 2D simulation results for the flow around a single cylinder (Re=200).
Data provenance
C
¯
D
C
D
′
C
l
′
St
Ding et al. [30]
1.348
0.035
0.466
0.196
Harichandan and Roy [25]
1.320
0.035
0.426
0.192
Lam et al. [17]
1.320
0.026
0.426
0.196
Mittal [31]
1.327
—
0.49
0.195
Current work
1.337
0.031
0.474
0.191Table 6
The numerical results at different Re.
Re
Data provenance
C
¯
D
C
l
′
St
300
SSTk-w model
1.327
0.558
0.203
Standardk-w model
1.341
0.554
0.200
Kravchenko and Moin: 3D [32]
1.28
0.4
0.203
Shen et al.: 2D [33]
1.33
—
0.21
Mittal: 2D [31]
1.357
0.652
0.208
500
Batcho and Karniadakis: 3D [34]
1.33
0.566
—
Tong et al.: 2D [20]
1.224
0.356
0.211
SSTk-w model
1.333
0.667
0.215
Standardk-w model
1.339
0.630
0.205
1000
Vu et al.: 2D [26]
1.259
—
0.2
Stringer et al.: 2D [23]
1.278
0.683
0.227
Jester and Kallinderis: 2D [35]
1.51
0.99
0.25
SSTk-w model
1.348
0.760
0.221
Standardk-w model
1.338
0.726
0.211In order to further verify the suitability of the calculation model and selected mesh for solving multicylinder problem, the comparison results for flow past two tandem cylinders atRe=200 and Re=1000 are presented. As shown in Table 7, the cases at L/D=2.0 and 4.0 are adopted to verify the accuracy. The upstream cylinder is marked as cylinder 1 and the downstream cylinder is marked as cylinder 2. Generally, the related calculated parameters are in good agreements with the results of Han et al. [18], Jester and Kallinderis [35], and Meneghini et al. [36], although there is a little disparity for the values of St2 at Re=1000 and L/D=2.0. The difference of numerical method and calculation mesh could cause the disparity. Generally, the difference for force coefficients obtained by different investigators is below 7.5%. For Re>200, to some extent, there is little disparity between 2D simulation and 3D simulation, especially for the average force coefficients. However, just as Jester and Kallinderis [35] pointed out the flow features between two kinds of simulation are nearly identical at Re=1000. Moreover, Stringer et al. [23] and Vu et al. [26] carried out a series of numerical studies to investigate flow pattern and flow dynamics for multicylinders based on the 2D simulations at Re≤1000 and obtained with satisfying results. Therefore, it is reasonable to speculate that the results of our 2D simulation at Re≤1000 could be adopted.Table 7
2D comparable simulation results for the flow around two tandem cylinders atRe=200 and Re=1000.
Data provenance
L
/
D
Cylinder 1
Cylinder 2
C
¯
D
1
C
l
1
′
S
t
1
C
¯
D
2
C
l
2
′
S
t
2
R
e
=
200
Han et al. [18]
2.0
1.041
—
0.132
-
0.199
—
0.132
Meneghini et al. [36]
2.0
1.03
—
0.13
-
0.17
—
0.13
Present
2.0
1.057
0.035
0.139
-
0.169
0.157
0.139
Meneghini et al. [36]
4.0
1.18
—
0.174
0.38
—
0.174
Present
4.0
1.226
0.527
0.173
0.388
1.112
0.173
R
e
=
1000
Vu et al. [26]
2.0
0.982
—
—
-
0.277
—
—
Jester and Kallinderis [35]
2.0
0.92
0.028
0.17
-
0.39
0.085
0.17
Present
2.0
0.927
0.020
0.196
-
0.374
0.072
0.196
Vu et al. [26]
4.0
1.288
—
—
0.321
—
—
Jester and Kallinderis [35]
4.0
1.45
—
—
0.17
—
—
Present
4.0
1.351
0.854
0.214
0.172
1.362
0.214


### 3.2. Flow Pattern Analysis

Figure4 shows the instantaneous vorticity contours for four square-arranged cylinders at different spacing ratios L/D=1.5, 2.0, 2.5, 3.5, 4.0, 5.0 and Reynolds numbers Re=200, 300, and 1000. At small spacing ratio L/D=1.5, as shown in Figures 4(a), 4(b), and 4(c), it can be seen that the downstream cylinders are completely engulfed by the outer shear layers separated from the upstream cylinders surface. This flow pattern is defined as the stable shielding flow pattern, which is almost consistent with the previous results obtained from Lam et al. [17], Han et al. [18], and Han et al. [37]. At such small spacing ratio, the adjacent effect is prominent for the fact that the vortex shedding pattern in wake flow is analogous to that of a single cylinder. However, with increasing of Reynolds number, the shielding length of the outer shear layer becomes shorter and the strength of gap flow between the upper and lower row cylinders is enhanced, resulting in complex interferences of the shear layers behind the downstream cylinders and smaller size of the vortex in the wake. In other words, the repulsion between the wake flows behind downstream cylinders is enhanced with increasing of Reynolds number. The gap flow has a significant effect on wake flow pattern, as Abbasi et al. [19] pointed out that the jets between the gaps strongly influence the wake interaction at different Reynolds numbers and gap spacing. When L/D increases to 2.0, as shown in Figure 4(e), the outer free shear layers separated from the upstream cylinders appear to be wiggled around downstream cylinders, defined as the wiggling shielding flow pattern [17, 18]. The outer free shear layers from the upstream cylinders reattach onto the downstream cylinders at L/D=2.5. Similarly, Han et al. [37] observed the typical wiggled shielding flow pattern behind the four circular cylinders at L/D=2.5 and Re=150. However, as shown in Figure 4(g), the vorticity contour plot indicates that the flow pattern at L/D=2.5 and Re=200 is on the stage of transition from wiggled shielding flow pattern to vortex shedding flow pattern, indicating that the variation of Reynolds number significantly influences the flow pattern characteristics. At the same time, the significant recirculation zone formed between cylinder 1 and cylinder 4 and the vortices shed from the downstream cylinders interfere and restrain the development of shear layers from upstream cylinders in Figures 4(g), 4(h), and 4(i). As a whole, the wake vortices behind cylinder 3 are prominent, and the width of wake flow behind cylinder 4 is narrower compared to that of cylinder 3, which could be attributed to the effect of biased gap flow. Comparing to the cases of L/D=2.5, it can be found that the positions of shed vortices in wake flow gradually approach the downstream cylinders with increasing of Reynolds number. With increasing of spacing ratio to 3.5, 4.0, and 5.0, the four-cylinder array could be approximately regarded as two rows of two tandem cylinders, and the surrounding shed vortices behind the upper and lower row cylinders are not affected by each other and move almost parallel in the streamwise direction. At L/D=4.0 and 5.0, the shear layers separated from the upstream cylinders are not reattached on the downstream cylinders and mature shed vortices can be observed behind each cylinder, which is consistent with previous results [17, 18]. In addition, the anti-phase vortex shedding pattern is observed at L/D=3.5, and the recirculation region between upstream and downstream cylinders gradually approaches the upstream cylinder with increasing of Reynolds number.Figure 4
Instantaneous vorticity distribution plots for 2D simulations atRe=200, 300, and 1000, and L/D=1.5, 2.0, 2.5, 3.5, 4.0, and 5.0.It is worthwhile to note that there is a slight phase difference between the vortices shed from the upper and lower cylinders atL/D=4.0; however, such phase difference appears to have little changes with increasing of Reynolds number. Meanwhile, the anti-phase vortex shedding at L/D=3.5 and 5.0 is still stable with increasing of Reynolds number. Therefore, it may be concluded that the phase difference is mainly related with spacing ratio L/D. Generally, the flow pattern at larger Reynolds number transforms faster than that at lower Reynolds number when the spacing ratio L/D is constant. For instance, when L/D=2.5, the flow patterns in Figures 4(g) and 4(h) are completely different at Re=200 and 300. As shown in Figure 4(h), the recirculation zone is close to the upstream cylinder, and the free shear layers separated from the upstream cylinders roll up, shed, and hardly reattach onto the downstream cylinders. However, for the similar case in Figure 4(g), the recirculation zone close to the downstream cylinders prevented the formation of vortices and forced the free shear layers from upstream cylinders to reattach onto downstream cylinders. Based on the results obtained from Figure 4, it can be concluded that the difference of the flow patterns at various spacing ratio is significant.In order to further investigate the effect of Reynolds number on the wake flow structures behind the four-cylinder array, instantaneous vorticity contours with variation of Reynolds number at constant spacing ratioL/D=3.0 are presented in Figure 5. At Reynolds number Re=100 (Figure 5(a)), typical wiggling shielding flow patterns are observed, and the inner free shear layers separated from upstream cylinders reattach onto downstream cylinder; however the outer free shear layers are not reattached onto the downstream cylinders but wiggled around the cylinders. When Re increases to 200, as shown in Figure 5(b), there is a large scale recirculation zone between the upstream and downstream cylinders, and the flow pattern appears to be still on the stage of transition from the wiggled shielding flow to the vortex shedding flow pattern. As presented in Figure 5(c), at Re=300, the free shear layers begin to roll up into relative mature vortices and it is interesting to note that the width of wake flow behind upper row cylinders is narrower than that behind the lower row cylinders, which can be also observed for two side-by-side arranged cylinders [2, 38, 39]. This special flow pattern is regarded as biased flow pattern. At Re=500 and 1000 (see Figures 5(d) and 5(e)), the anti-phase vortex shedding is observed and the wake flow behind the upper and lower row cylinders is paralleled and symmetric about the centerline between upper and lower cylinders. In other words, the fluctuating lifts applied on the cylinder 1 and cylinder 2 or cylinder 3 and cylinder 4 oscillate synchronously with 180° out-of-phase difference. More details about anti-phase shedding vortex could be referenced from the visualization results of Lam et al. [17] at L/D=5.0 and Re=200 or the particle simulation results of Farrant et al. [16] at L/D=5.0 and Re=200. Biased flow patterns with wide wake flow in one side and narrow wake flow in the other side are common for flow past two side-by-side cylinders at intermediate spacing ratios. However, similar phenomenon for the four square-arranged cylinders [13] is sparser and the reasons for causing biased wake flow and bistable phenomenon are yet unclear. Indeed, the vortex shedding patterns at L/D=2.5 and in Figure 5(c) can be summarized as the biased flow pattern in this study.Figure 5
Instantaneous vorticity distribution plots atL/D=3.0 and Re=100, 200, 300, 500, and 1000.
(a)
(b)
(c)
(d)
(e)In conclusion, the spacing ratio has a significant effect on different flow patterns in the simulation of flow past four squared-arranged cylinders. Moreover, different flow patterns could be observed for different Reynolds numbers, even at a constant spacing ratio, especially for small spacing ratios (L/D≤2.5). Therefore, the combined effects of spacing ratio and Reynolds number must be comprehensively considered in order to get comprehensive and accurate details about the flow past four-cylinder array.


### 3.3. Force and Strouhal Numbers Analysis

In the following section, the mean and root-mean-square (rms) values of drag and lift coefficient, as well as the Strouhal numbers, are presented and discussed. The mean drag coefficients are simplified asC¯D1,C¯D2,C¯D3, and C¯D4, where the subscripts “1”, “2”, “3”, and “4” denote the upstream and downstream cylinders, respectively (refer to Figure 1). Similar definitions are used for the root-mean-square drag coefficients CD1′,CD2′,CD3′,CD4′, the mean lift coefficients Cl¯1,Cl¯2,Cl¯3,Cl¯4, and the root-mean-square lift coefficients Cl1′, Cl2′, Cl3′, Cl4′.Figure6 shows the mean drag coefficients on four cylinders with variation of spacing ratios and Reynolds numbers. As shown in Figure 6(a), it can be seen that at Re=100 the mean drag coefficients on the upper and lower cylinders (cylinders 1 and 2, cylinders 3 and 4) are almost the same at L/D≥2.0, and there is only slight disparity between the upper and lower cylinders at L/D<2.0. This phenomenon may suggest that the interferences between the inner free shear layers become weakened with the increasing of the gap between the upper and lower cylinders. At the same time, compared with the results of Lam et al. [17], C¯D1 and C¯D2 are relatively lower, while C¯D3 and C¯D4 are higher, which are consistent with the similar results at Re=200 (see Figure 6(b)). Such disparity could be attributed to the difference of computational mesh. For Re=100, the value of C¯D1 or C¯D2 decreases with increasing of L/D, reaches the minimum at around L/D=3.0, and then increases again. Analogous results could be found from Lam et al. [17], and the spacing ratio corresponding to the minimum value of C¯D1 or C¯D2 is 1.0 smaller than that of Lam et al. [17]. In other words, L/D corresponding to the step rising of C¯D3 or C¯D4 is 1.0 smaller than that of Lam et al. [17]. Similar conclusion could be obtained at Re=200 (see Figure 6(b)); however the difference of the spacing ratios between this study and Lam et al. [17] is decreased to 0.5. Such phenomenon may be attributed to the following viewer. Different computational mesh and solving scheme may cause different flow pattern at the same spacing ratio. Compared with Lam et al. [17] the flow pattern transformation in this study leads in the spacing ratio at a constant Reynolds number. That is reason for the occurrence of above phenomenon. It is noted that the mean drag coefficient at Re=300 is very close to that at Re=500, which could be attributed to the insignificant transformation of flow pattern. In general, C¯D1 and C¯D2 slightly decrease with L/D varying from 1.5 to 2.0 at the five different Reynolds numbers. This could be explained as the fact that, at small spacing ratio, the forces are mainly applied on the upstream cylinders because the downstream cylinders are completely shielded by free shear layers from upstream cylinders. However, at L/D≥2.0, the free shear layers reattach onto the downstream due to the weak strength of recirculation zone between upstream and downstream cylinder, so the forces almost completely applied on the upstream cylinders are shared by the downstream cylinders. As a result, C¯D1 and C¯D2 decrease while C¯D3 and C¯D4 increase. Similar phenomenon can be found at high Reynolds numbers in Wang et al. [12] and Lam and Fang [7], which obtained the consistent results at Re=8.0×103 and Re=12.8×103, respectively, as presented in Figure 6(d). However, the trough L/D corresponding to minimum C¯D1 or C¯D2 varies with Reynolds numbers. For instance, the trough spacing ratios L/D are 3.0 and 2.5, respectively, at Re=100 and Re=12.8×103. With increases of spacing ratio to 5.0, the free shear layers from upstream cylinders roll up into shed vortices, and C¯D1 and C¯D2 are very close to that of a single cylinder. Based on the above discussions, it can be seen that C¯D1 and C¯D2 decrease at first, then increase again, and finally approach the value of a single cylinder. Basically, the mean drag coefficients on upstream cylinders at high Reynolds numbers are slightly lower than those at low Reynolds numbers, which could be analogous to that of a single cylinder. However, some values of C¯D1 and C¯D2 at Re=300 are abnormally large, which may be attributed to the interference when flow state changes from laminar model to turbulent model. The downstream cylinders are simultaneously influenced by the flow wake from the upstream cylinders and the wake flow interferences in the lee of downstream cylinders; namely, the environment of the downstream cylinders is more severe than that of the upstream cylinders. For Re=100, C¯D3 and C¯D4 completely increase with spacing ratio L/D and similar phenomenon can be observed in Wang et al. (2013) and Lam and Fang [7]. At Reynolds numbers Re≥200, C¯D3 increased and then decreased with increasing of spacing ratio. The reattachment of free shear layers from upstream cylinders results in the increasing of C¯D3. However, when the flow pattern switches into the vortex shedding flow pattern, the shear layers generated from downstream cylinders are affected by the vortices shed from upstream cylinders, which could be the reason for the decreasing of C¯D3. It is also noted that the spacing ratios corresponding to the peak values of C¯D3 and C¯D4 are distinguishing at different Reynolds numbers. The variation of C¯D3 and C¯D4 indicates that the alternation of Reynolds numbers results in different flow patterns and the Reynolds number is a significant factor affecting the flow pattern. When L/D=3.0, the values of C¯D3 at Re=200, 300, and 500 are 3.55, 4.02, and 3.35 times larger than that at Re=100, which clearly indicates that the key mechanical parameters change significantly when flow pattern transformation occurs. Meanwhile, C¯D3 is 2.85 times larger than C¯D4 at L/D=3.0 and Re=300. Such great difference could be attributed to the special flow characteristics that the gap flow is biased with a wide wake flow and a narrow wake flow behind the downstream cylinders (see Figure 5(c)). Sumner [2] studied the biased gap flow for two side-by-side circular cylinders and pointed out that the vortex shedding frequency and drag force on the cylinder immersed in narrow wake flow are higher than those in wide wake flow. Similar conclusion is obtained for the flow past four-cylinder array. When the Reynolds number increases to 500, the difference between C¯D3 and C¯D4 decreases, but C¯D3 is nearly 1.69 times larger than C¯D4. It could be speculated that the flow pattern transformation causes the decrease of difference between C¯D3 and C¯D4, but the exiting difference at Re=500 indicates the asymmetry of the pressure distribution on the downstream cylinders. At Re=1000, such difference between C¯D3 and C¯D4 decreases largely, as presented in Figure 6(d). In addition, C¯D3 and C¯D4 increase with increasing of Reynolds number at spacing ratio L/D varying from 2.0 to 2.5. One noticeable feature in Figure 6(d) is that C¯D4 is negative at Re=1000 and L/D<2.0. Similarly, Abbasi et al. [19] found that C¯D3 (C¯D4) is negative at gap ratio g=1.0 and Re=130 (Re=110). Moreover, Lam and Fang [7] concluded that C¯D3 and C¯D4 are all negative at Re=12.8×103 and L/D=1.6. However, it is not the case for Re<1000. Lam et al. [17] and Abbasi et al. [19] explained that, at high Re, the free shear layers separated from cylinder 1, after shielding behind cylinder 4 immediately, induce a strong backflow that applies a negative drag on cylinder 4; however, at low Reynolds numbers, the wake vortex or the backflow zone is formed far away from the downstream cylinder and the strength of backflow is too weak to produce a negative drag on cylinder 4.Figure 6
The mean drag coefficients versusL/D at different Re: experimental studies (Lam and Fang (1995) [7]; Wang et al. (2013) [12]); simulation studies (Lam et al. (2008) [17]); and present results: (a) C¯D of four cylinders at Re=100, (b) C¯D of four cylinders at Re=200, (c) C¯D of four cylinders at Re=300 and 500, and (d) C¯D of four cylinders at Re=1000.
(a)
(b)
(c)
(d)In general, the symmetry characteristic for the mean lift coefficients is clearly presented for upper and lower row cylinders in Figure7. In other words, Cl¯1 and Cl¯2 are equal in magnitude and opposite in symbol, as well as Cl¯3 and Cl¯4. Cl¯1 and Cl¯2 gradually tend to zero, with increasing of spacing ratio. This may indicate that the interferences between cylinder 1 and cylinder 2 are weakened with increasing of L/D and the free shear layers from upstream cylinders roll up and shed at large spacing ratios, which is similar to that of a single cylinder. As shown in Figures 7(a) and 7(b), it can be seen that the mean lift coefficients on the four cylinders at Re=100 and 200 are consistent with the result of Lam et al. [17]. By comparing Cl¯1 and Cl¯2 at different Reynolds numbers, it can be obtained that, at L/D≥2.5, the effects of Reynolds numbers on the magnitude of Cl¯1 or Cl¯2 are insignificant while at L/D<2.5 this kind of effect is significant and the smaller L/D is, the more significant such effect is. For example, when L/D is equal to 1.5, Cl¯1 at Re=100 is nearly 1.66 times larger than that at Re=1000. Moreover, the mean lift coefficient Cl¯1 of Lam et al. [17] at Re=200 is 4.9 times larger than that of Lam and Fang [7] at Re=12.8×103. These phenomena can also be associated with corresponding flow patterns and it can be concluded that at small spacing ratios (L/D<2.5) the mean lift coefficient is sensitive to the flow transformation. When the spacing ratio L/D reaches to 4.0, the effect of Reynolds numbers on Cl¯1 and Cl¯2 can be almost negligible. The vortex shedding behaviors of upstream cylinders are similar to that of an isolated circular cylinder, where Cl¯1 and Cl¯2 are almost zero.Figure 7
The mean lift coefficients versusL/D at different Re: experimental studies (Lam and Fang (1995) [7]; Wang et al. (2013) [12]); simulation studies (Lam et al. (2008) [17]); and present results: (a) C¯L of four cylinders at Re=100, (b) C¯L of four cylinders at Re=200, (c) C¯L of four cylinders at Re=300, (d) C¯L of four cylinders at Re=500, and (e) C¯L of four cylinders at Re=1000.
(a)
(b)
(c)
(d)
(e)Basically,Cl¯3 and Cl¯4 gradually tend to be zero with increasing of L/D and the magnitudes of Cl¯3 and Cl¯4 are obviously lower than those of Cl¯1 and Cl¯2. However, concerning the fluctuating characteristics of Cl¯3 and Cl¯4, it could be obtained that the lifts imposed on the downstream cylinders are sensitive to the flow pattern transformation. For instance, as shown in Figure 7(d), it is easy to find obvious concave and tortuous characteristics for curve of Cl¯3 varying with L/D at Re=500. The absolute values of Cl¯3 at L/D=2.0 and 3.0 are, respectively, 8.97 and 7.83 times larger than that at L/D=2.5, which may indicate the great difference among the three diverse flow patterns. Similarly, the concave characteristic can be found at Re=1000 with one interesting feature that Cl¯3 (Cl¯4) at L/D=2.0 and 3.0 are equal in magnitude and opposite in sign. However, reason for the interesting feature is unclear. Moreover, at Re=100, the absolute values of Cl¯3 and Cl¯4 appear to have no change when L/D≤3.0, which can be explained as the fact that the downstream cylinders are engulfed by the free shear layers separated from the upstream cylinders and the flow patterns are still on the stage of stable shielding flow pattern. However, at Re=200, Cl¯3 and Cl¯4 fluctuate significantly due to the obvious flow pattern transformation. In general, the mean lift coefficients of downstream cylinders are much lower than those of upstream cylinders regardless of the Reynolds number, since the downstream cylinders are always engulfed by the wake flow behind the upstream cylinders.Figure8 presents the nondimensional vortex shedding frequency Strouhal number for the four square-arranged cylinders at different Reynolds numbers and spacing ratios. The four primary Strouhal numbers are simplified as St1, St2, St3, and St4, where the subscripts “1”, “2”, “3”, and “4” denote the upstream and downstream cylinders, respectively (refer to Figure 1). It is difficult to figure out the primary vortex shedding frequency of each cylinder through the analysis of FFT (fast Fourier transform) for some cases, like Re=200 when L/D=1.5, so these corresponding Strouhal numbers are not presented in Figure 8. As a whole, the Strouhal numbers on four cylinders are almost the same at a constant spacing ratio and Reynolds number, which indicates the coupling vortex shedding characteristic between upstream and downstream cylinders. Since the amplitude of Strouhal numbers on four cylinders is almost the same, St1 is only described here. Generally speaking, for the cases L/D≤2.0, St1 differs obviously for different Reynolds numbers at a constant spacing ratio, while St1 is not sensitive to the variation of Reynolds number (Re≥200) at L/D>2.0. There is an evident concave characteristic for each curve of St regardless of the case Re=200; however, the spacing ratios corresponding to the trough value of St differ largely. Even so, St1 strictly increases with increasing of spacing ratio (L/D≥2.5) at any constant Reynolds number ranging from 200 to 1000. It is noted that, at Re=100, the bottom point of curve for St is far away from those at other Reynolds numbers for the fact that the flow pattern transformation at Re=100 is slow. When the spacing ratio L/D reaches to 5.0, Strouhal numbers on four cylinders are close to 0.2 except for the case at Re=100 and the effect of Reynolds number on Strouhal numbers can be neglected at such large spacing ratio. At higher Reynolds number Re=2100, St4 as reported in Lam and Lo [6] is cited for reference with L/D varying from 1.28 to 5.96 (see Figure 8(d)). At L/D>3.0, Strouhal numbers in this study are relatively higher than their resolution results; however, at L/D≤3.0, the vortices shed from cylinder 4 in particularly high frequency that significantly exceeds all the vortex frequencies computed in this study. Particularly, at L/D=1.28, St4 is up to 0.367, which could be attributed to the occurrence of the strong biased gap flow generated at small gap. Based on the comparable results presented in Figure 8, it is concluded that regardless of the Reynolds number (except for Re=200) Strouhal number decreases with the variation of L/D and then increases until it is close to the value of single cylinder at corresponding Reynolds number. And Strouhal numbers are sensitive to the flow pattern transformation at L/D≤3.0.Figure 8
The Strouhal numbers versusL/D at different Re: experimental studies (Lam and Lo (1992) [6]); simulation studies (Han et al. (2013) [18]; Lam et al. (2008) [17]); and present results: (a) St1 and St2 at Re=100 and 200, respectively, (b) St1 and St2 at Re=300, 500, and 1000, respectively, (c) St3 and St4 at Re=100 and 200, respectively, and (d) St3 and St4 at Re=300, 500, and 1000, respectively.
(a)
(b)
(c)
(d)The root-mean-square drag and lift coefficients, namely,CD′ and Cl′, are shown in Figures 9 and 10, respectively. The trend and amplitude of CD1′ and CD2′ (CD3′ and CD4′) are approximately the same. Similar conclusion can be obtained on Cl1′ and Cl2′ (Cl3′ and Cl4′). Significant increase or decrease of root-mean-square drag and lift coefficients indicated the occurrence of obvious flow pattern transformation. For instance, at Re=200, CD1′ or CD2′ increases at first, then decreases, and increases again at last. At the same time, Cl′ on four cylinders increase significantly at Re=200 and L/D ranging from 2.0 to 3.5. Particularly, at Re=1000, as shown in Figures 9(e) and 10(e), CD′ and Cl′ on four cylinders change significantly with spacing ratio increasing from 1.5 to 2.5. For instance, for the case of Re=1000, CD′ and Cl′ at L/D=2.0 are 3.15 and 4.08 times larger than those at L/D=1.5. Considering the instantaneous vorticity contours discussed above, it is easy to find that the server fluctuating characteristic of CD′ and Cl′ in various cases is closely related to the flow pattern transformation. Significant disparity for CD′ and Cl′ between this study and Lam et al. [17] at Re=100 and 200 could be found, especially at L/D≥3.0. At such critical spacing ratio (L/D≥3.0) the amplitude of the forces on cylinders may be sensitive to different computational mesh and solving scheme, which cause significant disparity for CD′ and Cl′ between this study and Lam et al. [17]. Basically, CD′ and Cl′ on four cylinders increase with increasing of Reynolds numbers from 100 to 1000 and appear to be larger than those reported in Lam et al. [8] and Wang et al. [12] at Re=4.1×104 and 8000. Perhaps, the flow features at such low Reynolds numbers (Re≤1000) are significantly different from the cases at relative high subcritical Reynolds numbers. It is also noted that CD′ and Cl′ on downstream cylinders are much larger than those on upstream cylinders due to the flow interferences for downstream cylinders. When the spacing ratio is larger than 4.0, CD′ and Cl′ on upstream cylinders gradually become stable, and the vortex shedding characteristic is similar to that of a single cylinder. CD3′ (Cl3′) and CD4′ (Cl4′) tend to be close to each other when L/D≥4.0. In Figure 9, significant spiking values of CD3′ and CD4′ can be found at various Reynolds ranging from 200 to 1000, which indicate the occurrence of significant flow pattern transformation. In other words, for these cases L/D= 2.0–3.0, CD′ on downstream cylinders sharply increase and then decrease at Re= 200–1000. However, the spacing ratios corresponding to the spiking values are different for different Reynolds numbers and the higher Reynolds number is, the lower spacing ratio L/D is. Similar phenomenon is also observed for Cl′ on downstream cylinders. Therefore, it can be obtained that the increase of Reynolds number can accelerate the flow pattern transformation.Figure 9
The root-mean-square drag coefficients versusL/D at different Re: experimental studies (Lam et al. (2003) [8]; Wang et al. (2013) [12]); simulation studies (Lam et al. (2008) [17]); and present results: (a) CD′ of four cylinders at Re=100, (b) CD′ of four cylinders at Re=200, (c) CD′ of four cylinders at Re=300, (d) CD′ of four cylinders at Re=500, and (e) CD′ of four cylinders at Re=1000.
(a)
(b)
(c)
(d)
(e)Figure 10
The root-mean-square lift coefficients versusL/D at different Re: experimental studies (Lam et al. (2003) [8]; Wang et al. (2013) [12]); simulation studies (Lam et al. (2008) [17]); and present results: (a) Cl′ of four cylinders at Re=100, (b) Cl′ of four cylinders at Re=200, (c)Cl′ of four cylinders at Re=300, (d) Cl′ of four cylinders at Re=500, and (e) Cl′ of four cylinders at Re=1000.
(a)
(b)
(c)
(d)
(e)


## 3.1. Validation Studies

In order to validate the accuracy of the present numerical model, two-dimensional simulation on flow past a single circular cylinder is carried out at different Reynolds numbersRe=200, 300, 500, and 1000. Figure 3 shows the schematic of the computational structured mesh and the local mesh near the cylinder surface. Four different structured finite element meshes based on the number of nodes on the cylinder surface are compared in Table 4. It can be seen that when the number of nodes on cylinder surface exceeds 140, the values of mean drag coefficient C¯D, root-mean-square drag, and lift coefficient CD′, Cl′ and Strouhal number St are almost stable. To reduce consumption of computational resources, the nodal point of 140 instead of 180 on the circumference of the cylinder surface was adopted in this study. The O-Grid technique is applied to refine the quality of mesh around the cylinder, where the minimum mesh size is 0.002D mm. The nondimensional time step (Δt=Ut/D) is 0.001 in all of the computational cases.Table 4
Grid independence test atRe=200.
Number of nodes on cylinder surface
Number of cells
Δ
t
C
¯
D
C
D
′
C
l
′
St
Grid 1
100
11328
0.001
1.330
0.030
0.466
0.193
Grid 2
140①
16048
0.001
1.337
0.031
0.474
0.191
Grid 3
180
20768
0.001
1.340
0.032
0.485
0.191
Grid 4
220
25488
0.001
1.340
0.032
0.486
0.192
①The final chosen calculation scheme.Figure 3
The view of structured mesh around the single circular cylinder. (a) Full view of the computation domain; (b) detail mesh close to the cylinder surface.
(a)
(b)Table5 shows the comparison of results for flow past a single cylinder from other numerical simulations at Re=200. For this Re, the present result of C¯D,CD′,Cl′, and St is consistent with the previous results obtained from Ding et al. [30], Harichandan and Roy [25], Lam et al. [17], and Mittal [31]. Compared with the result of Ding et al. [30], the value deviation is 0.82% for C¯D, 12.90% for CD′, 1.93% for Cl′, and 2.62% for St, respectively. When Reynolds number Re≥300, the comparison of results from different turbulent models is carried out. Previous numerical results for flow past single circular cylinder at different Re ranging from 300 to 1000 are presented in Table 6. Different simulation codes and grids may cause a bit disparity between other works [20, 23, 26, 31–35] for related parameters listed in Table 6. However, it can be seen from Table 6 that the results obtained from the Standardk-w model appear to be more consistent with the other numerical studies for the cases of Reynolds number Re=300, 500, and 1000.Table 5
Comparable 2D simulation results for the flow around a single cylinder (Re=200).
Data provenance
C
¯
D
C
D
′
C
l
′
St
Ding et al. [30]
1.348
0.035
0.466
0.196
Harichandan and Roy [25]
1.320
0.035
0.426
0.192
Lam et al. [17]
1.320
0.026
0.426
0.196
Mittal [31]
1.327
—
0.49
0.195
Current work
1.337
0.031
0.474
0.191Table 6
The numerical results at different Re.
Re
Data provenance
C
¯
D
C
l
′
St
300
SSTk-w model
1.327
0.558
0.203
Standardk-w model
1.341
0.554
0.200
Kravchenko and Moin: 3D [32]
1.28
0.4
0.203
Shen et al.: 2D [33]
1.33
—
0.21
Mittal: 2D [31]
1.357
0.652
0.208
500
Batcho and Karniadakis: 3D [34]
1.33
0.566
—
Tong et al.: 2D [20]
1.224
0.356
0.211
SSTk-w model
1.333
0.667
0.215
Standardk-w model
1.339
0.630
0.205
1000
Vu et al.: 2D [26]
1.259
—
0.2
Stringer et al.: 2D [23]
1.278
0.683
0.227
Jester and Kallinderis: 2D [35]
1.51
0.99
0.25
SSTk-w model
1.348
0.760
0.221
Standardk-w model
1.338
0.726
0.211In order to further verify the suitability of the calculation model and selected mesh for solving multicylinder problem, the comparison results for flow past two tandem cylinders atRe=200 and Re=1000 are presented. As shown in Table 7, the cases at L/D=2.0 and 4.0 are adopted to verify the accuracy. The upstream cylinder is marked as cylinder 1 and the downstream cylinder is marked as cylinder 2. Generally, the related calculated parameters are in good agreements with the results of Han et al. [18], Jester and Kallinderis [35], and Meneghini et al. [36], although there is a little disparity for the values of St2 at Re=1000 and L/D=2.0. The difference of numerical method and calculation mesh could cause the disparity. Generally, the difference for force coefficients obtained by different investigators is below 7.5%. For Re>200, to some extent, there is little disparity between 2D simulation and 3D simulation, especially for the average force coefficients. However, just as Jester and Kallinderis [35] pointed out the flow features between two kinds of simulation are nearly identical at Re=1000. Moreover, Stringer et al. [23] and Vu et al. [26] carried out a series of numerical studies to investigate flow pattern and flow dynamics for multicylinders based on the 2D simulations at Re≤1000 and obtained with satisfying results. Therefore, it is reasonable to speculate that the results of our 2D simulation at Re≤1000 could be adopted.Table 7
2D comparable simulation results for the flow around two tandem cylinders atRe=200 and Re=1000.
Data provenance
L
/
D
Cylinder 1
Cylinder 2
C
¯
D
1
C
l
1
′
S
t
1
C
¯
D
2
C
l
2
′
S
t
2
R
e
=
200
Han et al. [18]
2.0
1.041
—
0.132
-
0.199
—
0.132
Meneghini et al. [36]
2.0
1.03
—
0.13
-
0.17
—
0.13
Present
2.0
1.057
0.035
0.139
-
0.169
0.157
0.139
Meneghini et al. [36]
4.0
1.18
—
0.174
0.38
—
0.174
Present
4.0
1.226
0.527
0.173
0.388
1.112
0.173
R
e
=
1000
Vu et al. [26]
2.0
0.982
—
—
-
0.277
—
—
Jester and Kallinderis [35]
2.0
0.92
0.028
0.17
-
0.39
0.085
0.17
Present
2.0
0.927
0.020
0.196
-
0.374
0.072
0.196
Vu et al. [26]
4.0
1.288
—
—
0.321
—
—
Jester and Kallinderis [35]
4.0
1.45
—
—
0.17
—
—
Present
4.0
1.351
0.854
0.214
0.172
1.362
0.214


## 3.2. Flow Pattern Analysis

Figure4 shows the instantaneous vorticity contours for four square-arranged cylinders at different spacing ratios L/D=1.5, 2.0, 2.5, 3.5, 4.0, 5.0 and Reynolds numbers Re=200, 300, and 1000. At small spacing ratio L/D=1.5, as shown in Figures 4(a), 4(b), and 4(c), it can be seen that the downstream cylinders are completely engulfed by the outer shear layers separated from the upstream cylinders surface. This flow pattern is defined as the stable shielding flow pattern, which is almost consistent with the previous results obtained from Lam et al. [17], Han et al. [18], and Han et al. [37]. At such small spacing ratio, the adjacent effect is prominent for the fact that the vortex shedding pattern in wake flow is analogous to that of a single cylinder. However, with increasing of Reynolds number, the shielding length of the outer shear layer becomes shorter and the strength of gap flow between the upper and lower row cylinders is enhanced, resulting in complex interferences of the shear layers behind the downstream cylinders and smaller size of the vortex in the wake. In other words, the repulsion between the wake flows behind downstream cylinders is enhanced with increasing of Reynolds number. The gap flow has a significant effect on wake flow pattern, as Abbasi et al. [19] pointed out that the jets between the gaps strongly influence the wake interaction at different Reynolds numbers and gap spacing. When L/D increases to 2.0, as shown in Figure 4(e), the outer free shear layers separated from the upstream cylinders appear to be wiggled around downstream cylinders, defined as the wiggling shielding flow pattern [17, 18]. The outer free shear layers from the upstream cylinders reattach onto the downstream cylinders at L/D=2.5. Similarly, Han et al. [37] observed the typical wiggled shielding flow pattern behind the four circular cylinders at L/D=2.5 and Re=150. However, as shown in Figure 4(g), the vorticity contour plot indicates that the flow pattern at L/D=2.5 and Re=200 is on the stage of transition from wiggled shielding flow pattern to vortex shedding flow pattern, indicating that the variation of Reynolds number significantly influences the flow pattern characteristics. At the same time, the significant recirculation zone formed between cylinder 1 and cylinder 4 and the vortices shed from the downstream cylinders interfere and restrain the development of shear layers from upstream cylinders in Figures 4(g), 4(h), and 4(i). As a whole, the wake vortices behind cylinder 3 are prominent, and the width of wake flow behind cylinder 4 is narrower compared to that of cylinder 3, which could be attributed to the effect of biased gap flow. Comparing to the cases of L/D=2.5, it can be found that the positions of shed vortices in wake flow gradually approach the downstream cylinders with increasing of Reynolds number. With increasing of spacing ratio to 3.5, 4.0, and 5.0, the four-cylinder array could be approximately regarded as two rows of two tandem cylinders, and the surrounding shed vortices behind the upper and lower row cylinders are not affected by each other and move almost parallel in the streamwise direction. At L/D=4.0 and 5.0, the shear layers separated from the upstream cylinders are not reattached on the downstream cylinders and mature shed vortices can be observed behind each cylinder, which is consistent with previous results [17, 18]. In addition, the anti-phase vortex shedding pattern is observed at L/D=3.5, and the recirculation region between upstream and downstream cylinders gradually approaches the upstream cylinder with increasing of Reynolds number.Figure 4
Instantaneous vorticity distribution plots for 2D simulations atRe=200, 300, and 1000, and L/D=1.5, 2.0, 2.5, 3.5, 4.0, and 5.0.It is worthwhile to note that there is a slight phase difference between the vortices shed from the upper and lower cylinders atL/D=4.0; however, such phase difference appears to have little changes with increasing of Reynolds number. Meanwhile, the anti-phase vortex shedding at L/D=3.5 and 5.0 is still stable with increasing of Reynolds number. Therefore, it may be concluded that the phase difference is mainly related with spacing ratio L/D. Generally, the flow pattern at larger Reynolds number transforms faster than that at lower Reynolds number when the spacing ratio L/D is constant. For instance, when L/D=2.5, the flow patterns in Figures 4(g) and 4(h) are completely different at Re=200 and 300. As shown in Figure 4(h), the recirculation zone is close to the upstream cylinder, and the free shear layers separated from the upstream cylinders roll up, shed, and hardly reattach onto the downstream cylinders. However, for the similar case in Figure 4(g), the recirculation zone close to the downstream cylinders prevented the formation of vortices and forced the free shear layers from upstream cylinders to reattach onto downstream cylinders. Based on the results obtained from Figure 4, it can be concluded that the difference of the flow patterns at various spacing ratio is significant.In order to further investigate the effect of Reynolds number on the wake flow structures behind the four-cylinder array, instantaneous vorticity contours with variation of Reynolds number at constant spacing ratioL/D=3.0 are presented in Figure 5. At Reynolds number Re=100 (Figure 5(a)), typical wiggling shielding flow patterns are observed, and the inner free shear layers separated from upstream cylinders reattach onto downstream cylinder; however the outer free shear layers are not reattached onto the downstream cylinders but wiggled around the cylinders. When Re increases to 200, as shown in Figure 5(b), there is a large scale recirculation zone between the upstream and downstream cylinders, and the flow pattern appears to be still on the stage of transition from the wiggled shielding flow to the vortex shedding flow pattern. As presented in Figure 5(c), at Re=300, the free shear layers begin to roll up into relative mature vortices and it is interesting to note that the width of wake flow behind upper row cylinders is narrower than that behind the lower row cylinders, which can be also observed for two side-by-side arranged cylinders [2, 38, 39]. This special flow pattern is regarded as biased flow pattern. At Re=500 and 1000 (see Figures 5(d) and 5(e)), the anti-phase vortex shedding is observed and the wake flow behind the upper and lower row cylinders is paralleled and symmetric about the centerline between upper and lower cylinders. In other words, the fluctuating lifts applied on the cylinder 1 and cylinder 2 or cylinder 3 and cylinder 4 oscillate synchronously with 180° out-of-phase difference. More details about anti-phase shedding vortex could be referenced from the visualization results of Lam et al. [17] at L/D=5.0 and Re=200 or the particle simulation results of Farrant et al. [16] at L/D=5.0 and Re=200. Biased flow patterns with wide wake flow in one side and narrow wake flow in the other side are common for flow past two side-by-side cylinders at intermediate spacing ratios. However, similar phenomenon for the four square-arranged cylinders [13] is sparser and the reasons for causing biased wake flow and bistable phenomenon are yet unclear. Indeed, the vortex shedding patterns at L/D=2.5 and in Figure 5(c) can be summarized as the biased flow pattern in this study.Figure 5
Instantaneous vorticity distribution plots atL/D=3.0 and Re=100, 200, 300, 500, and 1000.
(a)
(b)
(c)
(d)
(e)In conclusion, the spacing ratio has a significant effect on different flow patterns in the simulation of flow past four squared-arranged cylinders. Moreover, different flow patterns could be observed for different Reynolds numbers, even at a constant spacing ratio, especially for small spacing ratios (L/D≤2.5). Therefore, the combined effects of spacing ratio and Reynolds number must be comprehensively considered in order to get comprehensive and accurate details about the flow past four-cylinder array.


## 3.3. Force and Strouhal Numbers Analysis

In the following section, the mean and root-mean-square (rms) values of drag and lift coefficient, as well as the Strouhal numbers, are presented and discussed. The mean drag coefficients are simplified asC¯D1,C¯D2,C¯D3, and C¯D4, where the subscripts “1”, “2”, “3”, and “4” denote the upstream and downstream cylinders, respectively (refer to Figure 1). Similar definitions are used for the root-mean-square drag coefficients CD1′,CD2′,CD3′,CD4′, the mean lift coefficients Cl¯1,Cl¯2,Cl¯3,Cl¯4, and the root-mean-square lift coefficients Cl1′, Cl2′, Cl3′, Cl4′.Figure6 shows the mean drag coefficients on four cylinders with variation of spacing ratios and Reynolds numbers. As shown in Figure 6(a), it can be seen that at Re=100 the mean drag coefficients on the upper and lower cylinders (cylinders 1 and 2, cylinders 3 and 4) are almost the same at L/D≥2.0, and there is only slight disparity between the upper and lower cylinders at L/D<2.0. This phenomenon may suggest that the interferences between the inner free shear layers become weakened with the increasing of the gap between the upper and lower cylinders. At the same time, compared with the results of Lam et al. [17], C¯D1 and C¯D2 are relatively lower, while C¯D3 and C¯D4 are higher, which are consistent with the similar results at Re=200 (see Figure 6(b)). Such disparity could be attributed to the difference of computational mesh. For Re=100, the value of C¯D1 or C¯D2 decreases with increasing of L/D, reaches the minimum at around L/D=3.0, and then increases again. Analogous results could be found from Lam et al. [17], and the spacing ratio corresponding to the minimum value of C¯D1 or C¯D2 is 1.0 smaller than that of Lam et al. [17]. In other words, L/D corresponding to the step rising of C¯D3 or C¯D4 is 1.0 smaller than that of Lam et al. [17]. Similar conclusion could be obtained at Re=200 (see Figure 6(b)); however the difference of the spacing ratios between this study and Lam et al. [17] is decreased to 0.5. Such phenomenon may be attributed to the following viewer. Different computational mesh and solving scheme may cause different flow pattern at the same spacing ratio. Compared with Lam et al. [17] the flow pattern transformation in this study leads in the spacing ratio at a constant Reynolds number. That is reason for the occurrence of above phenomenon. It is noted that the mean drag coefficient at Re=300 is very close to that at Re=500, which could be attributed to the insignificant transformation of flow pattern. In general, C¯D1 and C¯D2 slightly decrease with L/D varying from 1.5 to 2.0 at the five different Reynolds numbers. This could be explained as the fact that, at small spacing ratio, the forces are mainly applied on the upstream cylinders because the downstream cylinders are completely shielded by free shear layers from upstream cylinders. However, at L/D≥2.0, the free shear layers reattach onto the downstream due to the weak strength of recirculation zone between upstream and downstream cylinder, so the forces almost completely applied on the upstream cylinders are shared by the downstream cylinders. As a result, C¯D1 and C¯D2 decrease while C¯D3 and C¯D4 increase. Similar phenomenon can be found at high Reynolds numbers in Wang et al. [12] and Lam and Fang [7], which obtained the consistent results at Re=8.0×103 and Re=12.8×103, respectively, as presented in Figure 6(d). However, the trough L/D corresponding to minimum C¯D1 or C¯D2 varies with Reynolds numbers. For instance, the trough spacing ratios L/D are 3.0 and 2.5, respectively, at Re=100 and Re=12.8×103. With increases of spacing ratio to 5.0, the free shear layers from upstream cylinders roll up into shed vortices, and C¯D1 and C¯D2 are very close to that of a single cylinder. Based on the above discussions, it can be seen that C¯D1 and C¯D2 decrease at first, then increase again, and finally approach the value of a single cylinder. Basically, the mean drag coefficients on upstream cylinders at high Reynolds numbers are slightly lower than those at low Reynolds numbers, which could be analogous to that of a single cylinder. However, some values of C¯D1 and C¯D2 at Re=300 are abnormally large, which may be attributed to the interference when flow state changes from laminar model to turbulent model. The downstream cylinders are simultaneously influenced by the flow wake from the upstream cylinders and the wake flow interferences in the lee of downstream cylinders; namely, the environment of the downstream cylinders is more severe than that of the upstream cylinders. For Re=100, C¯D3 and C¯D4 completely increase with spacing ratio L/D and similar phenomenon can be observed in Wang et al. (2013) and Lam and Fang [7]. At Reynolds numbers Re≥200, C¯D3 increased and then decreased with increasing of spacing ratio. The reattachment of free shear layers from upstream cylinders results in the increasing of C¯D3. However, when the flow pattern switches into the vortex shedding flow pattern, the shear layers generated from downstream cylinders are affected by the vortices shed from upstream cylinders, which could be the reason for the decreasing of C¯D3. It is also noted that the spacing ratios corresponding to the peak values of C¯D3 and C¯D4 are distinguishing at different Reynolds numbers. The variation of C¯D3 and C¯D4 indicates that the alternation of Reynolds numbers results in different flow patterns and the Reynolds number is a significant factor affecting the flow pattern. When L/D=3.0, the values of C¯D3 at Re=200, 300, and 500 are 3.55, 4.02, and 3.35 times larger than that at Re=100, which clearly indicates that the key mechanical parameters change significantly when flow pattern transformation occurs. Meanwhile, C¯D3 is 2.85 times larger than C¯D4 at L/D=3.0 and Re=300. Such great difference could be attributed to the special flow characteristics that the gap flow is biased with a wide wake flow and a narrow wake flow behind the downstream cylinders (see Figure 5(c)). Sumner [2] studied the biased gap flow for two side-by-side circular cylinders and pointed out that the vortex shedding frequency and drag force on the cylinder immersed in narrow wake flow are higher than those in wide wake flow. Similar conclusion is obtained for the flow past four-cylinder array. When the Reynolds number increases to 500, the difference between C¯D3 and C¯D4 decreases, but C¯D3 is nearly 1.69 times larger than C¯D4. It could be speculated that the flow pattern transformation causes the decrease of difference between C¯D3 and C¯D4, but the exiting difference at Re=500 indicates the asymmetry of the pressure distribution on the downstream cylinders. At Re=1000, such difference between C¯D3 and C¯D4 decreases largely, as presented in Figure 6(d). In addition, C¯D3 and C¯D4 increase with increasing of Reynolds number at spacing ratio L/D varying from 2.0 to 2.5. One noticeable feature in Figure 6(d) is that C¯D4 is negative at Re=1000 and L/D<2.0. Similarly, Abbasi et al. [19] found that C¯D3 (C¯D4) is negative at gap ratio g=1.0 and Re=130 (Re=110). Moreover, Lam and Fang [7] concluded that C¯D3 and C¯D4 are all negative at Re=12.8×103 and L/D=1.6. However, it is not the case for Re<1000. Lam et al. [17] and Abbasi et al. [19] explained that, at high Re, the free shear layers separated from cylinder 1, after shielding behind cylinder 4 immediately, induce a strong backflow that applies a negative drag on cylinder 4; however, at low Reynolds numbers, the wake vortex or the backflow zone is formed far away from the downstream cylinder and the strength of backflow is too weak to produce a negative drag on cylinder 4.Figure 6
The mean drag coefficients versusL/D at different Re: experimental studies (Lam and Fang (1995) [7]; Wang et al. (2013) [12]); simulation studies (Lam et al. (2008) [17]); and present results: (a) C¯D of four cylinders at Re=100, (b) C¯D of four cylinders at Re=200, (c) C¯D of four cylinders at Re=300 and 500, and (d) C¯D of four cylinders at Re=1000.
(a)
(b)
(c)
(d)In general, the symmetry characteristic for the mean lift coefficients is clearly presented for upper and lower row cylinders in Figure7. In other words, Cl¯1 and Cl¯2 are equal in magnitude and opposite in symbol, as well as Cl¯3 and Cl¯4. Cl¯1 and Cl¯2 gradually tend to zero, with increasing of spacing ratio. This may indicate that the interferences between cylinder 1 and cylinder 2 are weakened with increasing of L/D and the free shear layers from upstream cylinders roll up and shed at large spacing ratios, which is similar to that of a single cylinder. As shown in Figures 7(a) and 7(b), it can be seen that the mean lift coefficients on the four cylinders at Re=100 and 200 are consistent with the result of Lam et al. [17]. By comparing Cl¯1 and Cl¯2 at different Reynolds numbers, it can be obtained that, at L/D≥2.5, the effects of Reynolds numbers on the magnitude of Cl¯1 or Cl¯2 are insignificant while at L/D<2.5 this kind of effect is significant and the smaller L/D is, the more significant such effect is. For example, when L/D is equal to 1.5, Cl¯1 at Re=100 is nearly 1.66 times larger than that at Re=1000. Moreover, the mean lift coefficient Cl¯1 of Lam et al. [17] at Re=200 is 4.9 times larger than that of Lam and Fang [7] at Re=12.8×103. These phenomena can also be associated with corresponding flow patterns and it can be concluded that at small spacing ratios (L/D<2.5) the mean lift coefficient is sensitive to the flow transformation. When the spacing ratio L/D reaches to 4.0, the effect of Reynolds numbers on Cl¯1 and Cl¯2 can be almost negligible. The vortex shedding behaviors of upstream cylinders are similar to that of an isolated circular cylinder, where Cl¯1 and Cl¯2 are almost zero.Figure 7
The mean lift coefficients versusL/D at different Re: experimental studies (Lam and Fang (1995) [7]; Wang et al. (2013) [12]); simulation studies (Lam et al. (2008) [17]); and present results: (a) C¯L of four cylinders at Re=100, (b) C¯L of four cylinders at Re=200, (c) C¯L of four cylinders at Re=300, (d) C¯L of four cylinders at Re=500, and (e) C¯L of four cylinders at Re=1000.
(a)
(b)
(c)
(d)
(e)Basically,Cl¯3 and Cl¯4 gradually tend to be zero with increasing of L/D and the magnitudes of Cl¯3 and Cl¯4 are obviously lower than those of Cl¯1 and Cl¯2. However, concerning the fluctuating characteristics of Cl¯3 and Cl¯4, it could be obtained that the lifts imposed on the downstream cylinders are sensitive to the flow pattern transformation. For instance, as shown in Figure 7(d), it is easy to find obvious concave and tortuous characteristics for curve of Cl¯3 varying with L/D at Re=500. The absolute values of Cl¯3 at L/D=2.0 and 3.0 are, respectively, 8.97 and 7.83 times larger than that at L/D=2.5, which may indicate the great difference among the three diverse flow patterns. Similarly, the concave characteristic can be found at Re=1000 with one interesting feature that Cl¯3 (Cl¯4) at L/D=2.0 and 3.0 are equal in magnitude and opposite in sign. However, reason for the interesting feature is unclear. Moreover, at Re=100, the absolute values of Cl¯3 and Cl¯4 appear to have no change when L/D≤3.0, which can be explained as the fact that the downstream cylinders are engulfed by the free shear layers separated from the upstream cylinders and the flow patterns are still on the stage of stable shielding flow pattern. However, at Re=200, Cl¯3 and Cl¯4 fluctuate significantly due to the obvious flow pattern transformation. In general, the mean lift coefficients of downstream cylinders are much lower than those of upstream cylinders regardless of the Reynolds number, since the downstream cylinders are always engulfed by the wake flow behind the upstream cylinders.Figure8 presents the nondimensional vortex shedding frequency Strouhal number for the four square-arranged cylinders at different Reynolds numbers and spacing ratios. The four primary Strouhal numbers are simplified as St1, St2, St3, and St4, where the subscripts “1”, “2”, “3”, and “4” denote the upstream and downstream cylinders, respectively (refer to Figure 1). It is difficult to figure out the primary vortex shedding frequency of each cylinder through the analysis of FFT (fast Fourier transform) for some cases, like Re=200 when L/D=1.5, so these corresponding Strouhal numbers are not presented in Figure 8. As a whole, the Strouhal numbers on four cylinders are almost the same at a constant spacing ratio and Reynolds number, which indicates the coupling vortex shedding characteristic between upstream and downstream cylinders. Since the amplitude of Strouhal numbers on four cylinders is almost the same, St1 is only described here. Generally speaking, for the cases L/D≤2.0, St1 differs obviously for different Reynolds numbers at a constant spacing ratio, while St1 is not sensitive to the variation of Reynolds number (Re≥200) at L/D>2.0. There is an evident concave characteristic for each curve of St regardless of the case Re=200; however, the spacing ratios corresponding to the trough value of St differ largely. Even so, St1 strictly increases with increasing of spacing ratio (L/D≥2.5) at any constant Reynolds number ranging from 200 to 1000. It is noted that, at Re=100, the bottom point of curve for St is far away from those at other Reynolds numbers for the fact that the flow pattern transformation at Re=100 is slow. When the spacing ratio L/D reaches to 5.0, Strouhal numbers on four cylinders are close to 0.2 except for the case at Re=100 and the effect of Reynolds number on Strouhal numbers can be neglected at such large spacing ratio. At higher Reynolds number Re=2100, St4 as reported in Lam and Lo [6] is cited for reference with L/D varying from 1.28 to 5.96 (see Figure 8(d)). At L/D>3.0, Strouhal numbers in this study are relatively higher than their resolution results; however, at L/D≤3.0, the vortices shed from cylinder 4 in particularly high frequency that significantly exceeds all the vortex frequencies computed in this study. Particularly, at L/D=1.28, St4 is up to 0.367, which could be attributed to the occurrence of the strong biased gap flow generated at small gap. Based on the comparable results presented in Figure 8, it is concluded that regardless of the Reynolds number (except for Re=200) Strouhal number decreases with the variation of L/D and then increases until it is close to the value of single cylinder at corresponding Reynolds number. And Strouhal numbers are sensitive to the flow pattern transformation at L/D≤3.0.Figure 8
The Strouhal numbers versusL/D at different Re: experimental studies (Lam and Lo (1992) [6]); simulation studies (Han et al. (2013) [18]; Lam et al. (2008) [17]); and present results: (a) St1 and St2 at Re=100 and 200, respectively, (b) St1 and St2 at Re=300, 500, and 1000, respectively, (c) St3 and St4 at Re=100 and 200, respectively, and (d) St3 and St4 at Re=300, 500, and 1000, respectively.
(a)
(b)
(c)
(d)The root-mean-square drag and lift coefficients, namely,CD′ and Cl′, are shown in Figures 9 and 10, respectively. The trend and amplitude of CD1′ and CD2′ (CD3′ and CD4′) are approximately the same. Similar conclusion can be obtained on Cl1′ and Cl2′ (Cl3′ and Cl4′). Significant increase or decrease of root-mean-square drag and lift coefficients indicated the occurrence of obvious flow pattern transformation. For instance, at Re=200, CD1′ or CD2′ increases at first, then decreases, and increases again at last. At the same time, Cl′ on four cylinders increase significantly at Re=200 and L/D ranging from 2.0 to 3.5. Particularly, at Re=1000, as shown in Figures 9(e) and 10(e), CD′ and Cl′ on four cylinders change significantly with spacing ratio increasing from 1.5 to 2.5. For instance, for the case of Re=1000, CD′ and Cl′ at L/D=2.0 are 3.15 and 4.08 times larger than those at L/D=1.5. Considering the instantaneous vorticity contours discussed above, it is easy to find that the server fluctuating characteristic of CD′ and Cl′ in various cases is closely related to the flow pattern transformation. Significant disparity for CD′ and Cl′ between this study and Lam et al. [17] at Re=100 and 200 could be found, especially at L/D≥3.0. At such critical spacing ratio (L/D≥3.0) the amplitude of the forces on cylinders may be sensitive to different computational mesh and solving scheme, which cause significant disparity for CD′ and Cl′ between this study and Lam et al. [17]. Basically, CD′ and Cl′ on four cylinders increase with increasing of Reynolds numbers from 100 to 1000 and appear to be larger than those reported in Lam et al. [8] and Wang et al. [12] at Re=4.1×104 and 8000. Perhaps, the flow features at such low Reynolds numbers (Re≤1000) are significantly different from the cases at relative high subcritical Reynolds numbers. It is also noted that CD′ and Cl′ on downstream cylinders are much larger than those on upstream cylinders due to the flow interferences for downstream cylinders. When the spacing ratio is larger than 4.0, CD′ and Cl′ on upstream cylinders gradually become stable, and the vortex shedding characteristic is similar to that of a single cylinder. CD3′ (Cl3′) and CD4′ (Cl4′) tend to be close to each other when L/D≥4.0. In Figure 9, significant spiking values of CD3′ and CD4′ can be found at various Reynolds ranging from 200 to 1000, which indicate the occurrence of significant flow pattern transformation. In other words, for these cases L/D= 2.0–3.0, CD′ on downstream cylinders sharply increase and then decrease at Re= 200–1000. However, the spacing ratios corresponding to the spiking values are different for different Reynolds numbers and the higher Reynolds number is, the lower spacing ratio L/D is. Similar phenomenon is also observed for Cl′ on downstream cylinders. Therefore, it can be obtained that the increase of Reynolds number can accelerate the flow pattern transformation.Figure 9
The root-mean-square drag coefficients versusL/D at different Re: experimental studies (Lam et al. (2003) [8]; Wang et al. (2013) [12]); simulation studies (Lam et al. (2008) [17]); and present results: (a) CD′ of four cylinders at Re=100, (b) CD′ of four cylinders at Re=200, (c) CD′ of four cylinders at Re=300, (d) CD′ of four cylinders at Re=500, and (e) CD′ of four cylinders at Re=1000.
(a)
(b)
(c)
(d)
(e)Figure 10
The root-mean-square lift coefficients versusL/D at different Re: experimental studies (Lam et al. (2003) [8]; Wang et al. (2013) [12]); simulation studies (Lam et al. (2008) [17]); and present results: (a) Cl′ of four cylinders at Re=100, (b) Cl′ of four cylinders at Re=200, (c)Cl′ of four cylinders at Re=300, (d) Cl′ of four cylinders at Re=500, and (e) Cl′ of four cylinders at Re=1000.
(a)
(b)
(c)
(d)
(e)


## 4. Conclusions

Numerical simulation on flow around four square-arranged cylinders is carried out to investigate the effect of Reynolds number and spacing ratio on wake flow structures. The flow characteristics such as vortex shedding pattern, force coefficients, and Strouhal number behind the four cylinders are discussed in this paper. Some conclusions are summarized as follows:(1)
Three basic flow patterns are observed with variation of spacing ratios, namely, the stable shielding flow pattern, wiggling shielding flow pattern, and vortex shedding flow pattern. The other flow patterns could be classified as the transition patterns between these basic flow patterns. The biased flow pattern with narrow and wide wakes is also observed atL/D=2.5 and L/D=3.0 (Re=300) for flow past four cylinders. It is found that the vortex shedding phase difference is only related with the spacing ratio and has little influence with Reynolds number. At L/D=1.5, the strength of gap flow is enhanced and the width of the vortices in wake decreases with increasing of Reynolds number. At L/D≤2.5, the flow interferences in the wake of downstream cylinders become more complex with increasing of Reynolds number.(2)
When the flow pattern switches, the values ofC¯D and Cl¯ on four cylinders change significantly. The spacing ratio L/D=3.0 is found to be closed to a critical spacing ratio for conspicuous flow pattern transformation at Re≥200 at which CD′ and Cl′ on four cylinders are increased significantly. Moreover, the fluctuating characteristic of the forces on the downstream cylinders is more evident than those on the upstream cylinders because of the more complex interferences between the downstream cylinders. At a constant spacing ratio, with increasing of Reynolds number, the amplitudes of CD′ and Cl′ increase due to the occurrence of flow pattern transitions.(3)
Regardless of the Reynolds numbers, atL/D≥1.5, the vortex shedding frequencies of the four cylinders are essentially the same. With variation of the spacing ratio, the Strouhal numbers decrease firstly and then increase until they are close to that of a single cylinder. However, when Re≥200, the variation of Strouhal number is found to be insensitive to Reynolds number at L/D>2.0.



---

*Source: 6381256-2017-03-27.xml*