Mathematical Problems in Engineering
Volume 2017 (2017), Article ID 6381256, 18 pages
https://doi.org/10.1155/2017/6381256

Research Article

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

Yang-yang Gao,¹ Chang-shan Yin,¹ Hao-qiang Zhang,² Kang Yang,¹ Xi-zeng Zhao,¹ and Zhilin Sun¹

¹Ocean College, Zhejiang University, Zhoushan 316021, China
²China Communications Construction Company, Beijing 100088, China

Correspondence should be addressed to Xi-zeng Zhao; xizengzhao@zju.edu.cn

Received 16 October 2016; Accepted 4 January 2017; Published 27 March 2017

Academic Editor: Nicolas Gourdain

Copyright © 2017 Yang-yang Gao et al. This is an open access article distributed under the Creative Commons Attribution License, which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.

Abstract.

In this paper, numerical simulations of flow past four square-arranged cylinders are carried out at different spacing ratios (; is the center to center distance; is the cylinder diameter) and Reynolds numbers . 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 .

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 cylinders ( is the center-center distance between cylinders, and 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 and (Reynolds number is defined as , is the diameter of cylinder, is the free-stream velocity, and is the kinematic viscosity) in an open-jet wind tunnel. The nondimensional vortex shedding frequency St (defined ; is the vortex shedding frequency) was found to be varied largely across wake at and equal to that for a single isolated cylinder when . Soon after that, Lam and Lo [6] categorized the flow characteristics behind four square-arranged cylinders into three flow patterns varied with ranging from 1.28 to 5.96 at : (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 ranges from 1.26 to 5.80 with and . They confirmed that is a critical spacing ratio and the forces characteristics vary largely with incidence angles (0°  ≤ α  ≤ 45°) and . Similarly, Lam et al. [8] experimentally measured the mean and fluctuating forces and Strouhal numbers on each cylinder at spacing ratio in the range 1.69 to 3.83, different angles of incidence α ranging from 0° to 180° at a 15° interval, and . As a result, three types of flow pattern were observed depending on : (i) shear layers around cylinder oscillate freely and vortex forms at large , (ii) shear layers reattach on or shield the downstream cylinder, and (iii) narrow gap flows form at small . Lam and Zou [9, 10] summarized several distinct flow patterns at subcritical Reynolds numbers 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 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 from 2.0 to 5.0, , 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 , , and (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 and the FEM study by Zhao and Cheng [15] at . Farrant et al. [16] employed the cell boundary element method (CBEM) to investigate the two-dimensional incompressible flow around four-cylinder array at 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 and 200. The results of their study showed that, depending on the spacing ratio , 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 with incidence angle = 0° and 45°. Three distinct flow patterns exist with variation of 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 () when spacing ratio was set at 2.0, 4.0, and 7.0, respectively. They also noted that the variation of was contributed much to the transformation of flow pattern and the characteristics of wake flow when spacing ratio is small (). Tong et al. [20] carried out a three-dimensional simulation on flow past four square-arranged cylinders at and . Four wake flow regimes were observed at different , 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

Technique

Experiments
Sayers [4, 5]
3.0 × 10⁴
1.1–5.0
0°–180°
Pressure method

Lam and Lo [6]
2.1 × 10³
1.28–5.96
0°–45°
LIF, PIV

Lam and Fang [7]
1.28 × 10⁴
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 × 10⁴
1.69–3.83
0°–180°
PELC

Lam and Zou [9]
(1.128–1.982) × 10⁴
1.5–5.0

LDA, PIV

Lam and Zou [10]
(1.1 –2.0) × 10⁴
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 × 10³–2.0 × 10⁴
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 × 10⁴
1.5, 3.5

LES

Lam et al. [17]
200
4.0

FVM

Zou et al. [14]
1.5 × 10⁴
1.5, 3.5

LES

Lam and Zou [10]
1.5 × 10⁴
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 ratio 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 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 () and subcritical flow () 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 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 for and 200, and the nondimensional tensorial forms of the Navier-Stokes governing equations for the incompressible viscous fluid flow can be expressed as

   are the velocities in the 2-dimensional Cartesian coordinates ( or ()), is the pressure, is the time, and is the Reynolds number. For the cases of , simulation could be conducted by directly solving (1).

When , 500, and 1000, the standard - 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

in (3) is Reynolds stress and it could be simplified as

when 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 and :

The definitions of related parameters in (3) and (5) are given as

with the corresponding explanation given as follows: , mean strain rate tensor; , Reynolds stress tensor. And the definitions of the other parameters are presented in the following:

, fluid density; , the molecular viscosity; , the Kronecker delta; , the eddy viscosity; , 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 is solved using (5), and (2) and (3) could be solved when the solution of 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 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:

2.2. Computational Domain and Boundary Conditions

Figure 1 shows the schematic of the geometry and mesh of flow past four square-arranged circular cylinders. The computational domain is ( is the diameter of cylinder) and the coordinate origin is set at the center position of the four-cylinder arrangement. The inlet of the computational domain is placed upstream from the coordinate origin and the outlet is placed downstream from the coordinate origin. The upper and lower boundary are placed 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
,  

Outflow
∂ ,  

Symmetry
V ,  

No-slip
,  

Table 3: Calculation models at different Reynolds numbers.

Re
(m/s)
Calculation model

100
0.005
Laminar

200
0.010
Laminar

300
0.015
Standard k-w

500
0.025
Standard k-w

1000
0.050
Standard k-w

(a)

(b)

Figure 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.

In this paper, laminar flow was simulated for the case of . At , 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 and St computed from laminar flow are close to the results computed from turbulent flow model [26]. The difference for 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 .

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 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⁻⁴ for continuity equation and 10⁻⁵ 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 numbers , 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 , root-mean-square drag, and lift coefficient , 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  mm. The nondimensional time step () is 0.001 in all of the computational cases.

Table 4: Grid independence test at .

Number of nodes on cylinder surface
Number of cells

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.

(a)

(b)

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.

Table 5 shows the comparison of results for flow past a single cylinder from other numerical simulations at . For this Re, the present result of , 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 , 12.90% for , 1.93% for , and 2.62% for St, respectively. When Reynolds number , the comparison of results from different turbulent models is carried out. Previous numerical results for flow past single circular cylinder at different 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 . However, it can be seen from Table 6 that the results obtained from the Standard k-w model appear to be more consistent with the other numerical studies for the cases of Reynolds number , 500, and 1000.

Table 5: Comparable 2D simulation results for the flow around a single cylinder ().

Data provenance

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.191

Table 6: The numerical results at different Re.

Re
Data provenance

St

300
SST k-w model
1.327
0.558
0.203

Standard k-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

SST k-w model
1.333
0.667
0.215

Standard k-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

SST k-w model
1.348
0.760
0.221

Standard k-w model
1.338
0.726
0.211

In 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 at and are presented. As shown in Table 7, the cases at 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 St₂ at and . 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 , 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 . 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 and obtained with satisfying results. Therefore, it is reasonable to speculate that the results of our 2D simulation at could be adopted.

Table 7: 2D comparable simulation results for the flow around two tandem cylinders at and .

Data provenance

Cylinder 1
Cylinder 2

Han et al. [18]
2.0
1.041
—
0.132

—
0.132

Meneghini et al. [36]
2.0
1.03
—
0.13

—
0.13

Present
2.0
1.057
0.035
0.139

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

Vu et al. [26]
2.0
0.982
—
—

—
—

Jester and Kallinderis [35]
2.0
0.92
0.028
0.17

0.085
0.17

Present
2.0
0.927
0.020
0.196

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

Figure 4 shows the instantaneous vorticity contours for four square-arranged cylinders at different spacing ratios , 2.0, 2.5, 3.5, 4.0, 5.0 and Reynolds numbers , 300, and 1000. At small spacing ratio , 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 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 . Similarly, Han et al. [37] observed the typical wiggled shielding flow pattern behind the four circular cylinders at and . However, as shown in Figure 4(g), the vorticity contour plot indicates that the flow pattern at and 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 , 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 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 , 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 at , 300, and 1000, and , 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 at ; however, such phase difference appears to have little changes with increasing of Reynolds number. Meanwhile, the anti-phase vortex shedding at 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 . Generally, the flow pattern at larger Reynolds number transforms faster than that at lower Reynolds number when the spacing ratio is constant. For instance, when , the flow patterns in Figures 4(g) and 4(h) are completely different at 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 ratio are presented in Figure 5. At Reynolds number (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 , 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 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 and or the particle simulation results of Farrant et al. [16] at and . 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 and in Figure 5(c) can be summarized as the biased flow pattern in this study.

(a)

(b)

(c)

(d)

(e)

Figure 5: Instantaneous vorticity distribution plots at and , 200, 300, 500, and 1000.

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 (). 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 as , and , 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 , the mean lift coefficients , and the root-mean-square lift coefficients , , , .

Figure 6 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 the mean drag coefficients on the upper and lower cylinders (cylinders 1 and 2, cylinders 3 and 4) are almost the same at , and there is only slight disparity between the upper and lower cylinders at . 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], and are relatively lower, while and are higher, which are consistent with the similar results at (see Figure 6(b)). Such disparity could be attributed to the difference of computational mesh. For , the value of or decreases with increasing of , reaches the minimum at around , and then increases again. Analogous results could be found from Lam et al. [17], and the spacing ratio corresponding to the minimum value of or is 1.0 smaller than that of Lam et al. [17]. In other words, corresponding to the step rising of or is 1.0 smaller than that of Lam et al. [17]. Similar conclusion could be obtained at (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 is very close to that at , which could be attributed to the insignificant transformation of flow pattern. In general, and slightly decrease with 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 , 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, and decrease while and 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 and , respectively, as presented in Figure 6(d). However, the trough corresponding to minimum or varies with Reynolds numbers. For instance, the trough spacing ratios are 3.0 and 2.5, respectively, at and . With increases of spacing ratio to 5.0, the free shear layers from upstream cylinders roll up into shed vortices, and and are very close to that of a single cylinder. Based on the above discussions, it can be seen that and 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 and at 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 , and completely increase with spacing ratio and similar phenomenon can be observed in Wang et al. (2013) and Lam and Fang [7]. At Reynolds numbers , increased and then decreased with increasing of spacing ratio. The reattachment of free shear layers from upstream cylinders results in the increasing of . 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 . It is also noted that the spacing ratios corresponding to the peak values of and are distinguishing at different Reynolds numbers. The variation of and 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 , the values of at , 300, and 500 are 3.55, 4.02, and 3.35 times larger than that at , which clearly indicates that the key mechanical parameters change significantly when flow pattern transformation occurs. Meanwhile, is 2.85 times larger than at and . 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 and decreases, but is nearly 1.69 times larger than . It could be speculated that the flow pattern transformation causes the decrease of difference between and , but the exiting difference at indicates the asymmetry of the pressure distribution on the downstream cylinders. At , such difference between and decreases largely, as presented in Figure 6(d). In addition, and increase with increasing of Reynolds number at spacing ratio varying from 2.0 to 2.5. One noticeable feature in Figure 6(d) is that is negative at and . Similarly, Abbasi et al. [19] found that () is negative at gap ratio and (). Moreover, Lam and Fang [7] concluded that and are all negative at and . However, it is not the case for . Lam et al. [17] and Abbasi et al. [19] explained that, at high , 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.

(a)

(b)

(c)

(d)

Figure 6: The mean drag coefficients versus at different : experimental studies (Lam and Fang (1995) [7]; Wang et al. (2013) [12]); simulation studies (Lam et al. (2008) [17]); and present results: (a) of four cylinders at , (b) of four cylinders at , (c) of four cylinders at and 500, and (d) of four cylinders at .

In general, the symmetry characteristic for the mean lift coefficients is clearly presented for upper and lower row cylinders in Figure 7. In other words, and are equal in magnitude and opposite in symbol, as well as and . and 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 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 and 200 are consistent with the result of Lam et al. [17]. By comparing and at different Reynolds numbers, it can be obtained that, at , the effects of Reynolds numbers on the magnitude of or are insignificant while at this kind of effect is significant and the smaller is, the more significant such effect is. For example, when is equal to 1.5, at is nearly 1.66 times larger than that at . Moreover, the mean lift coefficient of Lam et al. [17] at is 4.9 times larger than that of Lam and Fang [7] at . These phenomena can also be associated with corresponding flow patterns and it can be concluded that at small spacing ratios (<2.5) the mean lift coefficient is sensitive to the flow transformation. When the spacing ratio reaches to 4.0, the effect of Reynolds numbers on and can be almost negligible. The vortex shedding behaviors of upstream cylinders are similar to that of an isolated circular cylinder, where and are almost zero.

(a)

(b)

(c)

(d)

(e)

Figure 7: The mean lift coefficients versus at different : experimental studies (Lam and Fang (1995) [7]; Wang et al. (2013) [12]); simulation studies (Lam et al. (2008) [17]); and present results: (a) of four cylinders at , (b) of four cylinders at , (c) of four cylinders at , (d) of four cylinders at , and (e) of four cylinders at .

Basically, and gradually tend to be zero with increasing of and the magnitudes of and are obviously lower than those of and . However, concerning the fluctuating characteristics of and , 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 varying with at . The absolute values of at .0 and 3.0 are, respectively, 8.97 and 7.83 times larger than that at , which may indicate the great difference among the three diverse flow patterns. Similarly, the concave characteristic can be found at with one interesting feature that () at and 3.0 are equal in magnitude and opposite in sign. However, reason for the interesting feature is unclear. Moreover, at , the absolute values of and appear to have no change when , 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 , and 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.

Figure 8 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 St₁, St₂, St₃, and St₄, 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 when , 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, St₁ is only described here. Generally speaking, for the cases , St₁ differs obviously for different Reynolds numbers at a constant spacing ratio, while St₁ is not sensitive to the variation of Reynolds number () at . There is an evident concave characteristic for each curve of St regardless of the case ; however, the spacing ratios corresponding to the trough value of St differ largely. Even so, St₁ strictly increases with increasing of spacing ratio () at any constant Reynolds number ranging from 200 to 1000. It is noted that, at , 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 is slow. When the spacing ratio reaches to 5.0, Strouhal numbers on four cylinders are close to 0.2 except for the case at and the effect of Reynolds number on Strouhal numbers can be neglected at such large spacing ratio. At higher Reynolds number , St₄ as reported in Lam and Lo [6] is cited for reference with varying from 1.28 to 5.96 (see Figure 8(d)). At , Strouhal numbers in this study are relatively higher than their resolution results; however, at , the vortices shed from cylinder 4 in particularly high frequency that significantly exceeds all the vortex frequencies computed in this study. Particularly, at , St₄ 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 ) Strouhal number decreases with the variation of 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 .

(a)

(b)

(c)

(d)

Figure 8: The Strouhal numbers versus at different : experimental studies (Lam and Lo (1992) [6]); simulation studies (Han et al. (2013) [18]; Lam et al. (2008) [17]); and present results: (a) St₁ and St₂ at and 200, respectively, (b) St₁ and St₂ at , 500, and 1000, respectively, (c) St₃ and St₄ at and 200, respectively, and (d) St₃ and St₄ at , 500, and 1000, respectively.

The root-mean-square drag and lift coefficients, namely, and , are shown in Figures 9 and 10, respectively. The trend and amplitude of and ( and ) are approximately the same. Similar conclusion can be obtained on and ( and ). Significant increase or decrease of root-mean-square drag and lift coefficients indicated the occurrence of obvious flow pattern transformation. For instance, at , or increases at first, then decreases, and increases again at last. At the same time, on four cylinders increase significantly at and ranging from 2.0 to 3.5. Particularly, at , as shown in Figures 9(e) and 10(e), and on four cylinders change significantly with spacing ratio increasing from 1.5 to 2.5. For instance, for the case of , and at are 3.15 and 4.08 times larger than those at . Considering the instantaneous vorticity contours discussed above, it is easy to find that the server fluctuating characteristic of and in various cases is closely related to the flow pattern transformation. Significant disparity for and between this study and Lam et al. [17] at and 200 could be found, especially at . At such critical spacing ratio () the amplitude of the forces on cylinders may be sensitive to different computational mesh and solving scheme, which cause significant disparity for and between this study and Lam et al. [17]. Basically, and 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 and 8000. Perhaps, the flow features at such low Reynolds numbers () are significantly different from the cases at relative high subcritical Reynolds numbers. It is also noted that and 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, and on upstream cylinders gradually become stable, and the vortex shedding characteristic is similar to that of a single cylinder. () and () tend to be close to each other when . In Figure 9, significant spiking values of and 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 2.0–3.0, on downstream cylinders sharply increase and then decrease at 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 is. Similar phenomenon is also observed for on downstream cylinders. Therefore, it can be obtained that the increase of Reynolds number can accelerate the flow pattern transformation.

(a)

(b)

(c)

(d)

(e)

Figure 9: The root-mean-square drag coefficients versus at different : experimental studies (Lam et al. (2003) [8]; Wang et al. (2013) [12]); simulation studies (Lam et al. (2008) [17]); and present results: (a) of four cylinders at , (b) of four cylinders at , (c) of four cylinders at , (d) of four cylinders at , and (e) of four cylinders at .

(a)

(b)

(c)

(d)

(e)

Figure 10: The root-mean-square lift coefficients versus at different : experimental studies (Lam et al. (2003) [8]; Wang et al. (2013) [12]); simulation studies (Lam et al. (2008) [17]); and present results: (a) of four cylinders at , (b) of four cylinders at , (c) Cl ′ of four cylinders at , (d) of four cylinders at , and (e) of four cylinders at .

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 at and () 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 , the strength of gap flow is enhanced and the width of the vortices in wake decreases with increasing of Reynolds number. At , 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 of and on four cylinders change significantly. The spacing ratio is found to be closed to a critical spacing ratio for conspicuous flow pattern transformation at at which and 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 and increase due to the occurrence of flow pattern transitions.(3)Regardless of the Reynolds numbers, at , 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 , the variation of Strouhal number is found to be insensitive to Reynolds number at .

Competing Interests

The authors declare that there is no conflict of interests regarding the publication of this paper.

Acknowledgments

This research was supported by National Natural Science Foundation of China (Grant nos. 51409231, 51479175, and 51679212), Zhejiang Provincial Natural Science Foundation of China (Grant nos. LY14E090009 and LR16E090002), the Scientific Research Foundation for the Returned Overseas Chinese Scholars, State Education Ministry (no. 1685 ), and the State Key Laboratory of Ocean Engineering (Shanghai Jiao Tong University) (Grant no. 1312), China.

References

M. M. Zdravkovich, Flow around Circular Cylinders, University of Oxford, Oxford, UK, 2003.

D. Sumner, “Two circular cylinders in cross-flow: a review,” Journal of Fluids and Structures, vol. 26, no. 6, pp. 849–899, 2010.

P. W. Bearman, “Circular cylinder wakes and vortex-induced vibrations,” Journal of Fluids and Structures, vol. 27, no. 5-6, pp. 648–658, 2011.

A. T. Sayers, “Flow interference between four equispaced cylinders when subjected to a cross flow,” Journal of Wind Engineering and Industrial Aerodynamics, vol. 31, no. 1, pp. 9–28, 1988.

A. T. Sayers, “Vortex shedding from groups of three and four equispaced cylinders situated in a cross flow,” Journal of Wind Engineering and Industrial Aerodynamics, vol. 34, no. 2, pp. 213–221, 1990.

K. Lam and S. C. Lo, “A visualization study of cross-flow around four cylinders in a square configuration,” Journal of Fluids and Structures, vol. 6, no. 1, pp. 109–131, 1992.

K. Lam and X. Fang, “The effect of interference of four equispaced cylinders in cross flow on pressure and force coefficients,” Journal of Fluids and Structures, vol. 9, no. 2, pp. 195–214, 1995.

K. Lam, J. Y. Li, and R. M. C. So, “Force coefficients and Strouhal numbers of four cylinders in cross flow,” Journal of Fluids and Structures, vol. 18, no. 3-4, pp. 305–324, 2003.

K. Lam and L. Zou, “Experimental and numerical study for the cross-flow around four cylinders in an in-line square configuration,” Journal of Mechanical Science and Technology, vol. 21, no. 9, pp. 1338–1343, 2007.

K. Lam and L. Zou, “Experimental study and large eddy simulation for the turbulent flow around four cylinders in an in-line square configuration,” International Journal of Heat and Fluid Flow, vol. 30, no. 2, pp. 276–285, 2009.

K. Lam, J. Y. Li, K. T. Chan, and R. M. C. So, “Flow pattern and velocity field distribution of cross-flow around four cylinders in a square configuration at a low Reynolds number,” Journal of Fluids and Structures, vol. 17, no. 5, pp. 665–679, 2003.

X. K. Wang, K. Gong, H. Liu, J.-X. Zhang, and S. K. Tan, “Flow around four cylinders arranged in a square configuration,” Journal of Fluids and Structures, vol. 43, pp. 179–199, 2013.

L. Zou and Y. F. Lin, “Bistable flow Pattern for cylinder arrays with a small spacing ratio,” Journal of Visualization, vol. 12, no. 3, pp. 193–194, 2009.

L. Zou, Y.-F. Lin, and K. Lam, “Large-eddy simulation of flow around cylinder arrays at a subcritical reynolds number,” Journal of Hydrodynamics, vol. 20, no. 4, pp. 403–413, 2008.

M. Zhao and L. Cheng, “Numerical simulation of vortex-induced vibration of four circular cylinders in a square configuration,” Journal of Fluids and Structures, vol. 31, pp. 125–140, 2012.

T. Farrant, M. Tan, and W. G. Price, “A cell boundary element method applied to laminar vortex-shedding from arrays of cylinders in various arrangements,” Journal of Fluids and Structures, vol. 14, no. 3, pp. 375–402, 2000.

K. Lam, W. Q. Gong, and R. M. C. So, “Numerical simulation of cross-flow around four cylinders in an in-line square configuration,” Journal of Fluids and Structures, vol. 24, no. 1, pp. 34–57, 2008.

Z. Han, D. Zhou, X. Gui, and J. Tu, “Numerical study of flow past four square-arranged cylinders using spectral element method,” Computers and Fluids, vol. 84, pp. 100–112, 2013.

W. S. Abbasi, S. U. I. Islam, S. C. Saha, Y. T. Gu, and Z. C. Ying, “Effect of Reynolds numbers on flow past four square cylinders in an in-line square configuration for different gap spacings,” Journal of Mechanical Science and Technology, vol. 28, no. 2, pp. 539–552, 2014.

F. Tong, L. Cheng, M. Zhao, T. Zhou, and X.-B. Chen, “The vortex shedding around four circular cylinders in an in-line square configuration,” Physics of Fluids, vol. 26, no. 2, 2014.

D. C. Wilcox, “Reassessment of the scale-determining equation for advanced turbulence models,” AIAA Journal, vol. 26, no. 11, pp. 1299–1310, 1988.

J. Boussinesq, Essai sur la Théorie des Eaux Courantes, vol. 23 of Mémoires Présentées par Divers Savants à l'Académie des Sciences, Imprimerie Nationale, Paris, France, 1877.

R. M. Stringer, J. Zang, and A. J. Hillis, “Unsteady RANS computations of flow around a circular cylinder for a wide range of Reynolds numbers,” Ocean Engineering, vol. 87, pp. 1–9, 2014.

C.-N. Ji, W.-L. Chen, and W.-H. Xu, “Flow-induced vibrations of four square-arranged circular cylinders,” Journal of Vibration and Shock, vol. 35, no. 11, pp. 54–60, 2016 (Chinese).

A. B. Harichandan and A. Roy, “Numerical investigation of low Reynolds number flow past two and three circular cylinders using unstructured grid CFR scheme,” International Journal of Heat and Fluid Flow, vol. 31, no. 2, pp. 154–171, 2010.

H. C. Vu, J. Ahn, and J. H. Hwang, “Numerical simulation of flow past two circular cylinders in tandem and side-by-side arrangement at low Reynolds numbers,” KSCE Journal of Civil Engineering, vol. 20, no. 4, pp. 1594–1604, 2016.

W. Wang, W. Xu, Z. Zhai, Q. Meng, T. Chen, and Y. Luo, “Realization of PISO algorithm and comparison convergence characteristics with SIMPLE series algorithms,” Journal of North China Institute of Water Conservancy and Hydroelectric Power, vol. 28, no. 4, pp. 33–36, 2007 (Chinese).

W. Shyy, S. Thakur, and J. Wright, “Second-order upwind and central difference schemes for recirculating flow computation,” AIAA Journal, vol. 30, no. 4, pp. 923–932, 1992.

M.-J. Ni, W.-Q. Tao, and S.-J. Wang, “Stability-controllable second-order difference scheme for convection term,” Journal of Thermal Science, vol. 7, no. 2, pp. 119–130, 1998.

H. Ding, C. Shu, K. S. Yeo, and D. Xu, “Numerical simulation of flows around two circular cylinders by mesh-free least square-based finite difference methods,” International Journal for Numerical Methods in Fluids, vol. 53, no. 2, pp. 305–332, 2007.

S. Mittal, “Excitation of shear layer instability in flow past a cylinder at low Reynolds number,” International Journal for Numerical Methods in Fluids, vol. 49, no. 10, pp. 1147–1167, 2005.

A. Kravchenko and P. Moin, “B-spline methods and zonal grids for numerical simulations of turbulent flows,” Tech. Rep. TF-73, Flow Physics and Computational Division, Department of Mechanical Engineering, Stanford University, Stanford, Calif, USA, 1998.

L. Shen, M. Liu, L. Wu, P. Li, and D. Shu, “Numerical simulation of the flow around circular and square cylinders at sub-critical Reynolds number,” Journal of Waterway and Harbor, vol. 35, no. 3, pp. 228–229, 2014 (Chinese).

P. Batcho and G. E. Karniadakis, “Chaotic transport in two-and three-dimensional flow past a cylinder,” Physics of Fluids A, vol. 3, no. 5, pp. 1051–1061, 1991.

W. Jester and Y. Kallinderis, “Numerical study of incompressible flow about fixed cylinder pairs,” Journal of Fluids and Structures, vol. 17, no. 4, pp. 561–577, 2003.

J. R. Meneghini, F. Saltara, C. L. R. Siqueira, and J. A. Ferrari Jr., “Numerical simulation of flow interference between two circular cylinders in tandem and side-by-side arrangements,” Journal of Fluids and Structures, vol. 15, no. 2, pp. 327–350, 2001.

Z. Han, D. Zhou, Y. Chen, X. Gui, and J. Li, “Numerical simulation of cross-flow around multiple circular cylinders by spectral element method,” Acta Aerodynamica Sinica, vol. 32, no. 1, pp. 21–28, 2014 (Chinese).

S. Kumar, G. Laughlin, and C. Cantu, “Near-wake structure behind two circular cylinders in a side-by-side configuration with heat release,” Physical Review E, vol. 80, no. 6, Article ID 066307, 2009.

M. M. Alam and J. P. Meyer, “Two interacting cylinders in cross flow,” Physical Review E—Statistical, Nonlinear, and Soft Matter Physics, vol. 84, no. 5, Article ID 056304, 2011.