Numerical Study of Flow around Two Circular Cylinders in Tandem, Side-By-Side and Staggered Arrangements

Abstract

Simulations are presented for flow around pairs of circular cylinders at a Reynolds number of 3900. The flow is assumed to be two-dimensional and incompressible in nature and the simulations are performed using a RANS (Reynolds Averaged Navier Stokes) approach with a k-ε model. Simulations are performed for three different configurations of the cylinders: A tandem configuration where the line joining the centre of the cylinders is parallel to the mean flow direction; side-by-side, where the centre line is perpendicular to the mean flow direction; and staggered where the centre line is an angle α to the flow direction. Simulation results are presented for cylinder separations ranging from 1.125 to 4 diameters and for values of α between 10° and 60°. The results are presented and discussed in terms of the lift and drag coefficients, the Strouhal number, the vorticity field and the flow regimes observed. The results and flow regimes are also compared to previous observations at lower Reynolds numbers to investigate the Reynolds number dependence of the phenomena.

1. Introduction

The flow around a single circular cylinder has been investigated in detail and it is one of the classical problems of fluid mechanics with the dependence of the flow on the Reynolds number being well established [1,2]. However, when more than one cylinder is placed in the fluid flow, the entire behaviour of the flow becomes more complex, depending not only on the Reynolds number but also on the distance between the cylinders and their orientation to the flow direction. The flow around two circular cylinders has not been studied at the same depth as the flow around a single cylinder. For a two-cylinder configuration, the line joining the centre of the cylinders can be aligned with the flow direction in the following ways: Parallel to the flow direction giving a tandem arrangement; perpendicular to the flow direction giving a side-by-side orientation; or at an intermediate angle, giving a staggered arrangement.

Early studies of tandem configuration were performed by Igarashi [3] in the range of $8.7\times {10}^3\le Re\le 5.2\times {10}^4$and $1.03\le \frac{L}{D}\le 5.0$, where L is the centre-to-centre separation of the cylinders along the flow direction, D is the diameter and Re is the Reynolds number based on the diameter. $Re$ sensitivity was examined and six flow patterns were identified, which depended on $L$ and $Re$. These flow patterns were simplified to three regimes by Zhou and Yiu [4]: An “extended-body” regime at $1<\frac{L}{D}<2$ where the cylinders act as one bluff body; a “reattachment” regime at $2<\frac{L}{D}<5$ where instead of enclosing the downstream cylinder, the shear layers from the upstream cylinder reattach to the downstream cylinder; and a “co-shedding” regime at $\frac{L}{D}>5$, where $L$ is large enough for Kármán vortex streets to occur behind both cylinders. Sharman et al. [5] performed numerical predictions at $Re=100$ for $2D<L<10D$ and found $3.75<{\left(\frac{L}{D}\right)}_c<4$, where ${\left(\frac{L}{D}\right)}_c$ is the critical spacing between the cylinders, where a sudden significant change in the values of force coefficients and the Strouhal number occurred. Meneghini et al. [6] performed numerical simulation at $Re=200$ and found $3<{\left(\frac{L}{D}\right)}_c<4$.

In the side-by-side configuration, the $Re$ effect is less noticeable than in the tandem configuration, and three main flow patterns were identified [7]. In a “single-bluff-body” at $1<\frac{T}{D}<1.1-1.2$, where T is the centre-to-centre separation of the cylinders perpendicular to the flow direction, the cylinders act as a single bluff body. In “biased flow” at $1.1-1.2<\frac{T}{D}<2-2.2$, one cylinder has a narrow wake and higher values of shedding frequency, $f_s$, and drag coefficient, $C_d$, while the other cylinder has a wider wake and lower values of the same parameters. In “parallel vortex streets” at $\frac{T}{D}>2-2.2$, the flow field becomes symmetrical, and behind both cylinders, Kármán vortex streets occur with the same $f_s$. Ng and Ko [8] examined side-by-side cylinders numerically in the range of $100<Re<200$ for $\frac{T}{D}=1.25, 1.75, 3.0$ and $3.5$. Out-of-phase and in-phase vortex shedding and a maximum $C_d$ greater than in the single cylinder case were recorded. Sarvghad-Moghaddam et al. [9] used the finite volume approach at $Re=100, 200$ and ${10}^4$ for $1.5\le \frac{T}{D}\le 4$. For $Re={10}^4$, values of $\overline{C_d}$, the time-averaged drag coefficient, and the Strouhal number, $St$, for $\frac{T}{D}=3$ were found to be $1.24$ and $0.22$ and for $\frac{T}{D}=4$ were $1.12$ and $0.25$, respectively. Side-by-side configurations have also been considered recently at higher $Re$ experimentally by Alam et al. [10] ($Re=5.5\times {10}^4)$ and by Wu et al. [11] ($Re=1.4\times {10}^4)$using an LES simulation, with the focus of both studies being the flopping behaviour observed. Afgan et al. [12] also performed LES simulations on side-by-side cylinders, at $Re=3000$. Afgan et al. also observed flopping behaviour in the flow and at this Reynolds number and identified a “single-bluff-body” at $\frac{T}{D}~1$, “biased flow” at $1.25<\frac{T}{D}<1.75$, and “parallel vortex streets” at $\frac{T}{D}>2$. A 3D DNS simulation was presented by Chen et al. [13] at $Re=500$. Due to the 3D nature of this simulation, the effect of both the 3D structure of the vortices and the finite length of the cylinder could be assessed. Chen et al. determined that due to the transfer of energy from the initial 2D vortices into the third dimension, the vortices were not as strong as in the 2D case and the lift and drag coefficients and shedding frequencies were smaller, indicating that these could be overestimated in a 2D simulation. Additionally, Chen et al. [13] showed that for a cylinder greater than 6D in length, the end effects on the cylinder have minimal influence.

In the staggered configuration, Sumner et al. [14] used flow visualization and PIV in the range of $1.0\le \frac{P}{D}\le 5.0$, $0°\le \alpha \le 90°$ and $850<Re<190$, where $\alpha$ is the angle between the centre-to-centre line and the flow direction with $\alpha$ = 0 corresponding to a tandem orientation and P is the centre-to-centre separation. Nine flow patterns were identified and further grouped into three categories named “single bluff-body flow patterns”, “small-incidence-angle flow patterns” and “large-incidence-angle flow patterns”. Zhou et al. [15] classified four types of behaviour in the range of $1.2<\frac{P}{D}<4$ and $1.5\times {10}^3<Re<2\times {10}^4$. Type 1 occurred for $\frac{P}{D}<1.25$ and the cylinders act as a single body with $St$ independent of $Re$. In type 2, $\alpha <10°$, a single value of $St$was also observed, but the value was dependent on $Re$ with a sudden change occurring at a critical point in the range of $6.5\times {10}^3<Re<9.5\times {10}^3$, with the critical $Re$, which was dependent on $\frac{P}{D}$. Type 3 with $1.2-1.5<\frac{P}{D}<2.2$ and $10°<\alpha <75°$ had different values of $St$ for each cylinder, with the higher value of $St$ more sensitive to change in $Re$ than the other. In type 4, with $1.2-2.2<\frac{P}{D}$, the two cylinders have a different $St$ at low $Re$, where one is below the $St$ value for a single cylinder and the other is above. As Re increases, $St$for both cylinders gradually change towards the $St$ value of a single cylinder. Wong et al. [16] examined the reliance of flow classification on $Re$ in the range of $1.2<\frac{P}{D}<6$ for $Re=1.5\times {10}^3, 7\times {10}^3$ and $2\times {10}^4$ and found that the transition between different flow regimes have a significant dependency on $Re$. Akbari and Price [17] presented simulations of the staggered configuration at $Re=800$, established five flow patterns and obtained consistent results compared with existing experimental data.

Although numerical methods have gradually increased over the past two decades, the flow around two cylinders has not been fully explored due to its complexity. The majority of numerical methods in the literature are restricted to $Re\le 1000$, as opposed to a larger range in experimental studies [7,17,18]. Here we present simulation results for tandem, side-by-side and staggered configurations of two cylinders at $Re=3900$ within the subcritical regime where data exist for a single cylinder for validation purposes, but where there are limited data for two-cylinder configurations. Additionally, simulations are also presented at lower Reynold numbers to enable a comparison of two-cylinder configurations. The main aim of this work is to provide details of the flow regimes for different arrangements of two circular cylinders and to investigate the Reynolds number dependence of these regimes. This has applications in a wide range of air and water flows involving chimneys, skyscrapers, heat exchangers, bridges, piers, risers between offshore platforms and the ocean bed and electrical transmission lines.

Details of the numerical method applied in this study, including details of the computational grid, are set out in Section 2. In Section 3, the results are presented and discussed. Initially, simulation around a single cylinder is considered a validation case. The three orientations of tandem, side-by-side and staggered are then presented and discussed in Section 3.2.1, Section 3.2.2 and Section 3.2.3, respectively. Finally, Section 4 contains the conclusions drawn from the study. Additionally, a full set of the results for each configuration is presented in the Supplemental Material.

2. Numerical Method

Two-dimensional incompressible simulations were performed using CFD software STAR-CCM+, version 12.06.011. In the laminar regime, the fluid is described by the incompressible Continuity and Navier–Stokes equations:

$$\frac{\partial u_i}{\partial x_i}=0\mathrm{and}\frac{\partial u}{\partial t}+u_j\frac{\partial u_i}{\partial x_j}=-\frac{1}{\rho }\frac{\partial p}{\partial x_i}+\frac{\partial }{\partial x_j}\left(v\frac{\partial u_i}{\partial x_j}\right)$$

while in the turbulent regime, we applied the Reynolds Averaged Navier–Stokes (RANS) equations

$$\frac{\partial {\overline{u}}_i}{\partial x_i}=0\mathrm{and}\frac{\partial {\overline{u}}_i}{\partial t}+{\overline{u}}_j\frac{\partial {\overline{u}}_i}{\partial x_j}=-\frac{1}{\rho }\frac{\partial \overline{p}}{\partial x_i}+\frac{\partial }{\partial x_j}\left(v\frac{\partial {\overline{u}}_i}{\partial x_j}-\overline{{u^{\prime }}_i{u^{\prime }}_j }\right)$$

where $u$ and $p$ are the fluid velocity and pressure, $v$ is the kinematic viscosity, an overbar represents the time-averaged quantity, a prime represents the fluctuating component and summation is assumed over repeated indices. The Reynolds stresses, $\overline{{u^{\prime }}_i{u^{\prime }}_j }$ was modelled by the standard k-ε turbulent model, where

$$\overline{{u^{\prime }}_i{u^{\prime }}_j }=v_t\left(\frac{\partial {\overline{u}}_i}{\partial x_j}+\frac{\partial {\overline{u}}_j}{\partial x_i}\right)-\frac{2}{3}k{\delta }_{ij} ,$$

${\delta }_{ij}$ is the Kronecker delta, $k$ is the kinetic energy of the fluctuating motion and $v_t$ is the eddy-viscosity. This was modelled as

$$v_t=C_{\mu }\frac{k^2}{\epsilon }$$

where the kinetic energy and the dissipation rate, $\epsilon$, were modelled as

$$\frac{\partial k}{\partial t}+{\overline{u}}_i\frac{\partial k}{\partial x_i}=\frac{\partial }{\partial x_i}\left(\frac{{\nu }_t}{{\sigma }_k} \frac{\partial k}{\partial x_i}\right)+{\nu }_t \left(\frac{\partial {\overline{u}}_i}{\partial x_j}+\frac{\partial {\overline{u}}_j}{\partial x_i}\right)\frac{\partial {\overline{u}}_i}{\partial x_j}-\epsilon$$

$$\frac{\partial ϵ}{\partial t}+{\overline{u}}_i\frac{\partial \epsilon }{\partial x_i}=\frac{\partial }{\partial x_i}\left(\frac{{\nu }_t}{{\sigma }_{\epsilon }} \frac{\partial \epsilon }{\partial x_i}\right)+C_{1\epsilon } \frac{\epsilon }{k}{\nu }_t \left(\frac{\partial {\overline{u}}_i}{\partial x_j}+\frac{\partial {\overline{u}}_j}{\partial x_i}\right)\frac{\partial {\overline{u}}_i}{\partial x_j}-C_{2\epsilon } \frac{{\epsilon }^2}{k}$$

respectively, and the constants are given by $C_{\mu }=0.09$, $C_{1\epsilon }=1.44$, $C_{2\epsilon }=1.92$ and ${\sigma }_{\epsilon }=1.3$.

When simulating turbulent flow, the RANS approach adopted here is not the only option. Direct Numerical Simulation (DNS) [13] provides a direct solution of the Navier–Stokes equations without the need for a turbulent model, such as the k-ε turbulent model applied here. This is achieved by using a mesh, which is fine enough to capture the smallest length scales present in the turbulence. This comes at a computational cost, for example, the 3D simulations around a single cylinder at $Re=500$ presented in [13] took more than a month running on 576 cores with a speed of 2.0 GHz. Large Eddy Simulation (LES) [12] is an alternative approach, which requires an intermediate level of mesh density between that of RANS and DNS, enabling it to directly simulate more of the turbulent flow compared to RANS. LES can produce greater detail of the structures within the wave region, but this also comes at an additional computational cost with a single run for a single cylinder at $Re=3900$ reported to take $3.3\times {10}^5$ core hours [19]. For this reason, RANS is the standard approach in most industrial applications. The accuracy of the RANS approach depends on the turbulent model selected and the values used for the related constants. The performance of the turbulent models differs between applications, Reynolds number and even the required output, and the improvement and development of new turbulent models is an active area of research. In terms of flow past a cylinder, a number of authors have compared the different turbulent models [20,21]. Stringer et al. [20] suggest that the k-ε turbulent model can perform well at lower Reynolds numbers but suggest moving to a k-ω or SST k-ω model at higher Reynolds numbers. Rehman et al. [21] observed that, for the Reynold number considered here, the k-ε model accurately predicted the force coefficients but was less effective at visualising the vortex shedding.

In the present study, the rectangular fluid domain was created with dimensions of $50\mathrm{D}$ by $25\mathrm{D}$ where $\mathrm{D}$ is the diameter of the cylinder, which was 0.01 m.

A velocity inlet was applied at the upstream boundary with the velocity required to achieve the desired Reynold number. A pressure outlet was set at the downstream boundary and the other two boundaries were set as symmetry planes. The cylinder combination (or single cylinder) was placed with its centre 12D from the inlet and 12.5D from the symmetry planes and had no-slip boundary conditions.

A polygonal mesher was applied to the fluid domain with prism layers around the surface of the cylinder, and the base size of the final mesh was set to 0.0004 m. Four refined regions, VC1, VC2, VC3 and VC4, were set up with relative sizes of $50\%, 25\%, 12.5\%$ and $6.25\%$, respectively, compared to the base size, as shown in Figure 1. VC4 is the intersection of two circles of diameter 3D centred on the mid-points of the two cylinders. The dimensions of VC3 were determined such that the nearest point to a cylinder on the horizontal lines was 3.5D and on the vertical lines was 3D. The mesh consists of $\mathrm{110,733}$ cells, $\mathrm{329,967}$ interior faces and $\mathrm{219,978}$ vertices in the case of two cylinders, in the side-by-side configuration, where $\frac{T}{D}=2$. For other configurations, the same principle was applied to generate the mesh, but the number of cells, interior faces and vertices changes with the different arrangement of the two cylinders.

For the purpose of validating the accuracy of the results, a mesh convergence study was performed for the flow around a single cylinder at $Re=3900$, with respect to $\overline{C_d}$ and $St$. shown in

Table 1. Mesh convergence data for the flow around a cylinder at $Re=3900$.

Base Size [m]	$\overline{{\mathit{C}}_{\mathit{d}}}$	$\Delta \overline{{\mathit{C}}_{\mathit{d}}}$	$\mathit{S}\mathit{t}$	$\Delta \mathit{S}\mathit{t}$
$0.016$	$1.027$	$+7.878\%$	$0.217$	−0.05%
$0.008\times \sqrt{2}$	$1.026$	$+7.773\%$	$0.223$	−0.03%
$0.008$	$1.004$	$+5.462\%$	$0.225$	−0.02%
$0.004\times \sqrt{2}$	$0.996$	$+4.622\%$	$0.223$	−0.03%
$\mathbf{0.004}$	$\mathbf{0.958}$	$\mathbf{+}\mathbf{0.630}\mathbf{\%}$	$\mathbf{0.229}$	0.00%
$0.002\times \sqrt{2}$	$0.972$	$+2.101\%$	$0.229$	0.00%
$0.002$	$0.955$	$+0.315\%$	$0.227$	−0.01%
$0.001\times \sqrt{2}$	$0.952$	-	$0.229$	-

where $\Delta \overline{C_d}$ and $\Delta St$ represent the percentage differences compared to the finest mesh. Based on these results, a base size of 0.004 m was selected for further simulations. Full details of the mesh are shown in

Table 2. Details of the mesh controls and volume controls used to generate mesh.

Control	Value	Control	Value
Base Size	$0.004 \mathrm{m}$	Volumetric Control 1	$0.002 \mathrm{m}$
Surface Growth Rate	$1.06$	Volumetric Control 2	$0.001 \mathrm{m}$
Number of Prism Layers	$6$	Volumetric Control 3	$5.0\times {10}^{-4} \mathrm{m}$
Prism Layer Stretching	$1.2$	Volumetric Control 4	$2.5\times {10}^{-4} \mathrm{m}$
Prism Layer Total Thickness	$2.393\times {10}^{-4} \mathrm{m}$

and Figure 2.

Following Stringer [20], the prism layer parameters were selected to give Y+ < 1 and to expand to match the cell size in VC4. Figure 3 shows the variation in $Y+$ at $Re=3900$where the values are less than one everywhere, with a maximum value of 0.963. This configuration was applied for all $Re$ considered. The physics models used for the simulations at different $Re$ of the flow are presented in

Table 3. Physics models used for the simulations at different $\mathit{R}\mathit{e}$ of the flow.

Model	Reynolds Number
	2040	80,200	3900
Space	Two-Dimensional
Time	Steady	Implicit Unsteady
Material	Liquid
Flow	Segregated Flow
Equation of State	Constant Density
Viscous Regime	Laminar		Turbulent
RANS Model	N/A		k-ε
Near Wall Treatment	N/A		Enhanced Wall Treatment

.

The liquid used in the simulation was water with the constant density $\rho =1000 \frac{kg}{m^3}$ and dynamic viscosity $\mu =8.9\times {10}^{-4}Pa s$. The temporal discretization of the implicit unsteady solver is set to 2nd-order. The time step was determined based on

$$\mathsf{\Delta }t=0.2\frac{\pi D}{nU}$$

where n is the number of cells around the cylinder and U is the speed required to achieve the desired Reynolds number. This gives $\mathsf{\Delta }t=0.001 s$ for $Re=3900$. The turbulence specification used in the case of $Re=3900$ is k-ε, following Rahman et al. [21]. The inlet conditions for the turbulent model were $k=5.854$ × 10⁻⁴ $\frac{J}{kg}$ and $\epsilon$= $3.325$ × 10⁻⁴ $\frac{J}{kg. s}$

3. Results and Discussion

Simulations were performed for the configurations shown in

Table 4. Details of simulation configurations.

Case	$\mathit{R}\mathit{e}$	Type	L/T/P	α(deg)	Case	$\mathit{R}\mathit{e}$	Type	L/T/P	α(deg)
C1	20	Single	-	-	C14	200	S-by-S	T = 2D	-
C2	40	Single	-	-	C15	3900	S-by-S	T= 2D	-
C3	80	Single	-	-	C16	200	S-by-S	T = 3D	-
C4	200	Single	-	-	C17	3900	S-by-S	T = 3D	-
C5	3900	Single	-	-	C18	3900	Staggered	P = 1D	30
C6	200	Tandem	L = 1.5D	-	C19	3900	Staggered	P = 1D	60
C7	3900	Tandem	L = 1.5D	-	C20	3900	Staggered	P = 1.125D	60
C8	200	Tandem	L = 3D	-	C21	3900	Staggered	P = 1.5D	10
C9	3900	Tandem	L = 3D	-	C22	3900	Staggered	P = 1.5D	20
C10	200	Tandem	L = 4D	-	C23	3900	Staggered	P = 4D	10
C11	3900	Tandem	L = 4D	-	C24	3900	Staggered	P = 1.5D	30
C12	200	S-by-S	T = 1.5D	-	C25	3900	Staggered	P = 1.5D	60
C13	3900	S-by-S	T = 1.5D	-	C26	3900	Staggered	P = 2.5D	50

.

A full set of results in terms of contour maps of the vorticity field, drag and lift coefficient from each cylinder and, where the frequency is not regular, Fourier transforms, are included in the Supplemental Material. Here only the specific results that are pertinent to the discussion are included.

3.1. Single Cylinder

Prior to considering two-cylinder configurations, a range of simulations were run on a single cylinder to validate the numerical approach.

The results of simulations of flow around a single cylinder were as expected, based on the information available in the literature, with $C_d$ decreasing and $St$ increasing as $Re$ increases. Time averaging was applied to the fully developed drag coefficient due to observed fluctuations, the full details of which are presented in the Supplemental Material, where $\overline{C_d}$ is the time-averaged value $C_d$.

A comparison of $C_d$ and $St$ between the present study and literature data is presented in

Table 5. Comparison between experimental and numerical results from literature and present study for ${\mathit{C}}_{\mathit{d}}$ and $\mathit{S}\mathit{t}$ in the flow around single cylinder.

Case	C1	C2	C3		C4		C5
$Re$	$20$	$40$	$80$		$200$		3900
Parameters	$C_d$	$C_d$	$\overline{C_d}$	$St$	$\overline{C_d}$	$St$	$\overline{C_d}$	$St$
Exp. Tritton [25]	$2.22$	$1.65$	$1.32$	$0.154$	-	-	-	-
Ding et al. [26]	$2.14$	$1.58$	-	-	-	-	-	-
Franke et al. [27]	-	-	$1.35$	$0.152$	$1.31$	$0.194$	-	-
Braza et al. [28]	-	-	$1.39$	$0.200$	$1.39$	$0.200$	-	-
Meneghini et al. [6]	-	-	-	-	$1.30$	$0.196$	-	-
Exp. Norberg [22]	-	-	-	-	-	-	$0.98\pm 0.05$	$0.215\pm 0.005$
2D Rahman et al. [21]	-	-	-	-	-	-	$0.997$	$0.222$
3D LES Rajani et al. [23]	-	-	-	-	-	-	$1.01-1.05$	$0.210-0.214$
3D LES Ouvrard et al. [24]	-	-	-	-	-	-	$0.90-1.02$	$0.218-0.228$
Present Study	$2.13$	$1.58$	$1.40$	$0.157$	$1.36$	$0.198$	$0.958$	$0.229$

and shows good agreement, indicating the suitability of the model and mesh configurations used in the simulations. For $Re=3900$, the $\overline{C_d}$ in the present study is within the error range of the experimental data [22], and the results from the present study differ from the 2D RANS study by Rahman et al. [21] by 4%. It is approximately 5% smaller than the lowest value found in the 3D LES simulation by Rajani et al. [23] and within the range observed in 3D LES simulations by Ouvrard et al. [24]. The range of values for the LES simulations was due to different models. In terms of $St,$ its value is larger than the experimental range [22], 2D simulation [21] and 3D LES simulation [23,24] by 4%, 3% 7%, and 0.05%, respectively. This further suggests that the approximations introduced by performing the simulations in 2D and using a RANS approach are acceptable.

3.2. Dual Cylinders

All tandem and side-by-side arrangements were performed at $Re=3900$ and also at $Re=200$ for comparative purposes with numerical reference solutions $($Meneghini et al. 2001). The staggered arrangements are performed at $Re=3900$ and are grouped into categories for comparative purposes with the experimental work of Sumner et al. [14] (in the range of $850<Re<1900$).

3.2.1. Tandem Arrangement

Table 6. Comparison between numerical literature results at $Re=200$ and present study at $Re=200$ and $Re=3900$ for $C_d$ and $St$ for two cylinders in the tandem arrangement.

Gap L	1.5 D				3.0 D				4.0 D
Parameters	$\overline{C_{d_{UC}}}$	$\overline{C_{d_{DC}}}$	$S{t}_{UC}$	$S{t}_{DC}$	$\overline{C_{d_{UC}}}$	$\overline{C_{d_{DC}}}$	$S{t}_{UC}$	$S{t}_{DC}$	$\overline{C_{d_{UC}}}$	$\overline{C_{d_{DC}}}$	$S{t}_{UC}$	$S{t}_{DC}$
Meneghini et al. [6]	$1.06$	$-0.18$	$0.167$	$0.167$	$1.00$	$-0.08$	$0.125$	$0.125$	$1.18$	$0.38$	$0.174$	$0.174$
$\mathrm{Present}\mathrm{Study}\mathrm{at}Re=200$	$1.09$	$-0.20$	$0.170$	$0.170$	$1.02$	$-0.12$	$0.129$	$0.129$	$1.31$	$0.59$	$0.182$	$0.182$
$\mathrm{Present}\mathrm{Study}\mathrm{at}Re=3900$	$0.82$	$-0.28$	$0.203$	$0.203$	$0.93$	$0.54$	$0.191$	$0.191$	$0.93$	$0.37$	$0.204$	$0.204$

shows $\mathrm{St}$ and $\overline{C_d}$ for the Downstream Cylinder (DC) and Upstream Cylinder (UC).

In Figure 4, C6 and C7 represent the tandem configuration at $\frac{L}{D}=1.5$. C6 is in the “extended-body” regime, where two cylinders act as one bluff body. The DC is fully submerged in the wake of the UC. In C6, the separated shear layers from the UC are enclosing the DC, without any reattachment onto its surface before the shedding of vortices occurs in the wake of DC, and the gap between cylinders is considered to contain stagnant fluid [7] (Figure 4a). Whereas, in C7, the separated shear layers slightly reattach to the DC, and the fluid in the gap between cylinders is no longer stagnant (Figure 4b). The wake in both cases is narrower and eddies are less compacted than in the single-cylinder cases at the same $Re$. The $\overline{C_d}$ in both cases is positive for the UC and negative for the DC with lower values as $Re$ increases. The $C_d$ fluctuations increase with Re, but in both cases, the DC fluctuations are higher. Further simulations at lower values of $\frac{L}{D}$ at $Re=3900$ (see Supplemental Material) were found to remain in the “reattachment” regime, indicating that the “extended body” regime does not exist at this Reynolds number.

C8 and C9 represent the tandem configuration at $\frac{L}{D}=3$. C8 is in the “reattachment” regime, where two cylinders are positioned with enough distance between them that the separated shear layers from the UC cannot enclose the DC, but they reattach to the DC [7]. The reattachment of the separated shear layers from the UC to the DC is shown in Figure 5a.

Whereas in C9, vortex shedding is evident from the UC, which indicates the transition to the “co-shedding” regime, where two cylinders are sufficiently far apart for the vortex shedding to occur from UC and DC, as shown in Figure 5b. In C8, the wake of the DC still maintains the pattern observed in C6 at $\frac{L}{D}=1.5$. Whereas in C9, the flow pattern further downstream began to resemble the single-cylinder case. In C8, the $\overline{C_d}$ of the UC has decreased (compared to C6), but $\overline{C_d}$ for the DC had increased, while still maintaining a negative value. The negative value of $\overline{C_d}$ indicates that the DC is still within a low-pressure region created by the separated shear layers from the UC [6]. In C9, the $\overline{C_d}$ of UC is higher than at $\frac{L}{D}=1.5$, and $\overline{C_d}$ for the DC has become positive, which is in agreement with $\overline{C_d}$ approaching the single-cylinder case as the gap increases. The behaviour of the $C_d$ fluctuations is similar to the case of $\frac{L}{D}=1.5$, but in each cylinder, they are higher.

C10 and C11 represent a tandem configuration at $\frac{L}{D}=4$. Both cases are now in the transition to the “co-shedding” regime. The UC is shedding vortices in both cases, but the DC is not located entirely outside the UC vortex formation [7], as shown in Figure 6. The eddies in the wake of the DC in C11 are distorted and weaker in terms of magnitude in comparison to C9, which is one of the characteristics of further transition to the “co-shedding” regime [7].

In C10, the $\overline{C_d}$ of UC had increased to a value close to the single-cylinder case and the $\overline{C_d}$ of the DC is positive. Whereas in C11, the $\overline{C_d}$ of the UC maintains the same as at the lower $\frac{L}{D}$, and $\overline{C_d}$ of the DC has decreased, in contrast to the same arrangement at $Re=200$.

The results at $Re=200$ in C6, C8 and C10 are as expected in comparison to the work performed by Meneghini et al. [6]. The critical value of the gap, where the $C_d$ started to rise, was found to be in the range of $3<{\left(\frac{L}{D}\right)}_c<4$, which agrees with the data available in the literature [6]. The results of $Re=3900$ in C7, C9 and C11 demonstrate the sensitivity of the tandem arrangement to $Re$. The ${\left(\frac{L}{D}\right)}_c$ was found to be in the range of $1.5<{\left(\frac{L}{D}\right)}_c<3$. The simulation at $\frac{L}{D}=1.5$ is in the “reattachment” regime and simulations at $\frac{L}{D}=3$ and $\frac{L}{D}=4$ are in the transition to the “co-shedding” regime. Additional simulations at lower $\frac{L}{D}$ (see Supplemental Material) were unable to find a lower limit to the “reattachment” regime indicating that the “extended body” regime does not exist at $Re=3900$. This indicates that as Re increases, the range of the gap corresponding to the simplified flow patterns from the work performed by Zhou and Yiu [4] also increases. The $St$ in the tandem arrangement are the same for UC and DC in every simulation, and $\frac{L}{D}=3$ displayed the smallest values of $St$.

3.2.2. Side-by-Side Arrangement

The comparison of $C_d$ and $St$ between the present study and literature data in the case of a side-by-side arrangement is presented in

Table 7. Comparison between numerical literature results at $\mathit{R}\mathit{e}\mathbf{=}\mathbf{200}$ and $\mathit{R}\mathit{e}\mathbf{=}\mathbf{3900}$ and the present study at $Re=200$ and $Re=3900$ for $C_d$ and $St$ for two cylinders, in side-by-side arrangement.

Gap T	1.5 D				2.0 D				3.0 D
Parameters	$\overline{C_{d_{HC}}}$	$\overline{C_{d_{LC}}}$	$S{t}_{HC}$	$S{t}_{LC}$	$\overline{C_{d_{HC}}}$	$\overline{C_{d_{LC}}}$	$S{t}_{HC}$	$S{t}_{LC}$	$\overline{C_{d_{HC}}}$	$\overline{C_{d_{LC}}}$	$S{t}_{HC}$	$S{t}_{LC}$
Meneghini et al. [6] $Re=200$	1.32	1.32	-	-	1.42	1.42	-	-	1.34	1.34	$0.174$	$0.174$
$\mathrm{Present}\mathrm{Study}\mathrm{at}Re=200$	1.55	1.57	0.201	0.204	1.58	1.63	0.216	0.223	1.56	1.56	0.214	0.214
$\mathrm{Present}\mathrm{Study}\mathrm{at}Re=3900$	1.13	0.96	0.261	0.258	1.22	1.15	0.244	0.241	1.10	1.10	0.238	0.238
Afgan et al. [12]$Re=3000$	1.4	1.2	0.39	0.11	1.3	1.3	0.24	0.22	1.4	1.4	0.23	0.23

. Agreement with Meneghini et al. [6] is not as close as in the tandem case. The mesh density applied in this study is greater than that used in [6] by approximately an order of magnitude, suggesting that the mesh applied here was appropriate. The side-by-side simulations at $Re=200$ were re-run using a refined mesh in VC3 and VC4; however, the results did not change in a significant manner, confirming that the mesh is appropriate. We also note that in the tandem case, only a small discrepancy was observed, suggesting that the discrepancy here is not mesh-dependent. We note that the force data in this arrangement is less regular than in the tandem arrangement, which is discussed below and is evident in the Supplemental Material. It is possible that this may be the cause of the poorer agreements with [6], either due to differences in the averaging approach or a finer mesh being required in [6] to fully capture the phenomenon. We note that the $\overline{C_d}$ (and $\overline{C_l}$) obtained by Meneghini et al. [6] at $\frac{T}{D}=1.5$ and $\frac{T}{D}=2$ are the same for Higher Cylinder (HC) and Lower Cylinder (LC) in each case, which is in contrast to the present study and Sumner [7]. The $Re=3900$ results are also compared to Afgan et al. [12], where the data have been extracted from Figures 14 and 15 in [12]. For the Strouhal number, there is reasonable agreement at T = 2D and T = 3D; however, there is a significant discrepancy at T = 1.5D where Afgan et al. [12] reports a significantly different value for the two cylinders. There are differences between the mean drag coefficient reported here at $Re=3900$ and by Afgan et al. [12] at $Re=3000$; however, they are constant with the decrease in the coefficient with increasing Reynolds number observed in all the results in the table between $Re=200$ and $Re=3900$.

C12 and C13 represent a side-by-side configuration at $\frac{T}{D}=1.5$. Both cases are in the “biased flow” regime, where one cylinder has a narrower wake and higher values of $C_d$ and $St$compared to the other cylinder [7]. In C12, the narrower wake belongs to the LC, but it is not very noticeable from the vorticity contours shown in Figure 7a. Furthermore, the shedding of vortices occurs close to the cylinders, and frontal stagnation points have shifted towards the direction of the gap. In C13, the narrower wake belongs to the HC, which is clearly visible in Figure 7b, and the shedding of vortices is further downstream. In C12, the combined wake is characterized by strong vortices from the outer shear layers, and the deflection of the gap flow changes the direction of the vortices [7]. In C13, the combined wake is not deflected to the same extent, but eddies are distorted.

The ${\mathit{C}}_{\mathit{d}}$ and ${\mathit{C}}_{\mathit{l}}$ plots for both cases are shown in Figure 8 and Figure 9, respectively. They provide evidence of the “flopping” behaviour of the wake, which was described in experimental work performed by Kim and Durbin [29]. The wake is deflected towards one cylinder, which causes its drag to increase [6]. The difference in ${\mathit{C}}_{\mathit{d}}$ values between the HC and the LC in C12 is insignificant, as opposed to C13. Whereas in C13, the fluctuations of ${\mathit{C}}_{\mathit{d}}$ are lower and more predictable than in C12. In the configuration of $\frac{\mathit{T}}{\mathit{D}}=\mathbf{1.5}$, there is a repulsive force acting on both cylinders, which causes $\overline{{\mathit{C}}_{\mathit{l}}}$ to be negative for LC and positive for HC in both cases [6]. The amplitude of ${\mathit{C}}_{\mathit{l}}$ is low in C13, which can be explained by inspecting the vorticity contours in Figure 7b. The wake downstream of both cylinders is almost symmetrical before the shedding occurs further downstream.

Fourier analysis was performed in Matlab R2018a and is shown in Figure 10. The $\overline{f_s}$ value for C12 HC, C12 LC, C13 HC and C13 LC was found to be $0.3586 \mathrm{Hz}$, $0.3638 \mathrm{Hz}$, $9.0761\mathrm{Hz}$ and $8.9659\mathrm{Hz}$, respectively. The differences between HC and LC are not significant, but the general concept agrees with the work of Sumner [7].

C14 and C15 represent a side-by-side configuration at $\frac{T}{D}=2$. Both cases are still in the “biased flow” regime. The vorticity contours for C14 and C15 are shown in Figure 11a and Figure 11b, respectively. In C14, a smaller deflection of the frontal stagnation points towards the gap can be observed. The repulsive force is also smaller, but the wake in both cases is still combined and disorganised [6]. Shedding of vortices in C15 occurs close to the cylinder, unlike in C13.

The “flopping” behaviour of the wake is still present (see Supplemental Material) and shows evidence of a beat-like phenomenon. The values of $\overline{C_d}$ have increased for every cylinder, which agrees with Meneghini et al. [6].

C16 and C17 represent a side-by-side configuration at $\frac{T}{D}=3$. Both cases are in the “parallel vortex streets” regime, and the vorticity contours are shown in Figure 12. The “flopping” was not present for either cylinder in each case, and symmetry was restored, with each cylinder acting more as a single cylinder. However, the interference effects of close proximity still occur. The wakes of both cylinders in each case are synchronized in an anti-phase shedding pattern [7].

The $C_d$ fluctuations are identical for HC and LC in each case, along with $\overline{C_d}$. The solutions stabilised but did not regain fluctuations that occur in the single-cylinder case, in agreement with Meneghini et al. [6].

The side-by-side arrangement considered in this study at $Re=200$, in C12, C14 and C16, generally agrees with the work performed by Meneghini et al. [6]. The results at $Re=3900$ in C13, C15 and C17 show evidence that the side-by-side configuration is not as sensitive to changes in $Re$ as was found for the tandem configuration and are comparable with the work of Afgan et al. [12], which found the “biased flow” regime at $1.25\le \frac{T}{D}\le 1.75$ and the “parallel vortex streets” regime at $\frac{T}{D}\ge 2$ at $Re=3000$.

3.2.3. Staggered Arrangement

The values of $\overline{C_d}$ and $St$ for the staggered arrangements are presented in

Table 8. Numerical results of $C_d$ and $St$ for single bluff-body flow patterns from present study at $Re=3900$ in the flow around two cylinders, in staggered arrangement.

P, α, Case	$1.0 \mathit{D}, 30°,$C18				$1.0 \mathit{D}, 60°,$C19				$1.125 \mathit{D}, 60°,$C20
Parameter	$\overline{C_{d_{UC}}}$	$\overline{C_{d_{DC}}}$	$S{t}_{UC}$	$S{t}_{DC}$	$\overline{C_{d_{UC}}}$	$\overline{C_{d_{DC}}}$	$S{t}_{UC}$	$S{t}_{DC}$	$\overline{C_{d_{UC}}}$	$\overline{C_{d_{DC}}}$	$S{t}_{UC}$	$S{t}_{DC}$
Value	$0.88$	$0.69$	$0.146$	$0.146$	$1.57$	$1.72$	$0.117$	$0.117$	$1.42$	$0.96$	$0.123$	$0.123$
P, α, case	$1.5 \mathit{D}, 10°$, C21				$1.5 \mathit{D}, 20°,$C22				$4.0 \mathit{D}, 10°$, C23
Parameter	$\overline{C_{d_{UC}}}$	$\overline{C_{d_{DC}}}$	$S{t}_{UC}$	$S{t}_{DC}$	$\overline{C_{d_{UC}}}$	$\overline{C_{d_{DC}}}$	$S{t}_{UC}$	$S{t}_{DC}$	$\overline{C_{d_{UC}}}$	$\overline{C_{d_{DC}}}$	$S{t}_{UC}$	$S{t}_{DC}$
Value	$0.86$	$0.03$	$0.190$	$0.190$	$0.86$	$0.39$	$0.165$	$0.165$	$0.92$	$0.61$	$0.205$	$0.205$
P, α, case	$\mathbf{1.5} \mathit{D}\mathbf{,} \mathbf{30}\mathbf{°}$, C24				$\mathbf{1.5} \mathit{D}\mathbf{,} \mathbf{60}\mathbf{°}\mathbf{,}$C25				$\mathbf{2.5} \mathit{D}\mathbf{,} \mathbf{50}\mathbf{°}$, C26
Parameter	$\overline{C_{d_{UC}}}$	$\overline{C_{d_{DC}}}$	$S{t}_{UC}$	$S{t}_{DC}$	$\overline{C_{d_{UC}}}$	$\overline{C_{d_{DC}}}$	$S{t}_{UC}$	$S{t}_{DC}$	$\overline{C_{d_{UC}}}$	$\overline{C_{d_{DC}}}$	$S{t}_{UC}$	$S{t}_{DC}$
Value	$0.90$	$0.62$	$0.161$	$0.161$	$1.13$	$0.93$	$0.309$	$0.157$	$0.96$	$0.95$	$0.294$	$0.207$

. C18, C19 and C20 represent staggered configurations for “single bluff-body flow patterns”. C18 is a type 1 single bluff-body (SBB1), C19 is a type 2 single bluff-body (SBB2) and C20 is a base-bleed (BB). C18 is characterised by different lengths of the shear layers. The shear layers from the UC are significantly longer than the ones from the DC. The stretched shear layers are more likely to develop instabilities in the vortex shedding. In C19, the length of the shear layers is similar for both cylinders and the instabilities occur from the UC and the DC, simultaneously. In C20, there is a narrow gap between the cylinders, which allows fluid to flow through the base region. The gap causes stretching of the near-wake region, which causes the vortices to occur further downstream compared to Sumner et al. [14]. The vorticity contours of C18, C19 and C20 are shown in Figure 13a, Figure 13b and Figure 13c, respectively.

C21, C22 and C23 represent staggered configurations of “small-incidence-angle flow patterns”. C21 is a shear layer reattachment (SLR), C22 is an induced separation (IS) and C23 is a vortex impingement (VI). In C21, there is a gap between the cylinders, but the flow through that gap is prevented due to the separated shear layer from the upper side of the UC. This shear layer immediately reattaches onto the DC, which prevents the flow from penetrating the gap. As a result, the flow pattern of the combined wake resembles C18, where the lower shear layer of the UC stretches along the wake downstream, but with weaker magnitude in terms of vorticity and with smaller eddies. In C22, the DC moved further away from the wake of the UC, which allows the upper shear layer of the UC to be directed into the gap between the cylinders. This allows the flow to penetrate the gap and indicates the end of the stretched lower shear layer for the UC from the previous configuration. In C23, the DC is far enough from the UC to allow the shedding of vortices to occur behind both cylinders, which is similar to the “co-shedding” regime in the tandem configuration. The UC started to transition to the behaviour of the single-cylinder case, and the vortices shed from the UC impinge upon the DC [14]. The vorticity contours of C21, C22 and C23 are shown in Figure 14a, Figure 14b and Figure 14c, respectively.

In C23, the fluctuations for $C_d$ have two different peak-to-peak values that repeat every other period, as shown in Figure 15. This is different to the observations for the other cases (see Supplemental Material) and is a sign of instability most likely caused by the vortex shedding from both cylinders with the DC remaining in the wake of the UC. It is worth mentioning that this is not observed for the tandem configuration at $\frac{L}{D}=4$ at the same $Re$, suggesting that this phenomenon is dependent on the angle $\alpha$. No instability was observed in $C_l$.

C24, C25 and C26 represent the staggered configuration of “large-incidence-angle flow patterns”. C24 is a vortex pairing and enveloping pattern (VPE), C25 is a vortex pairing, splitting and enveloping pattern (VPSE) and C26 is a synchronized vortex shedding pattern (SVS). In C24, the inner shear layers of both cylinders and the vortices created from them further downstream are greatly synchronized, and this phenomenon is called the pairing of gap vortices, which are then enveloped by the vortices formed from the outer shear layer of the UC. In C25, the angle of incidence was doubled, which caused the enveloping to be unfinished due to the increased gap and therefore increased distance between the inner shear layers of both cylinders. The vortices created from those further separated inner shear layers are not well synchronized, as opposed to C24. Gap vortices created from each cylinder split into two different regions of vorticity, which disrupts the enveloping of the outer shear layer of the UC [14]. The vorticity contours of C24 and C25 are shown in Figure 16a and Figure 16b, respectively.

The $C_d$ results for C25 and, to a lesser extent, for C26 (see Supplemental Material) show significantly unstable and unpredictable fluctuations from the UC and the DC, providing evidence of “flopping” behaviour, which was already discussed in the side-by-side arrangement, at $\frac{T}{D}=1.5$ and $\frac{T}{D}=2$. Fourier analysis of both cylinders in C25 and C26 is shown in Figure 17. The $\overline{f_s}$ for C25 UC, C25 DC, C26 UC and C26 DC was found to be $10.7209 \mathrm{Hz}$, $5.4590 \mathrm{Hz}$, $10.2091\mathrm{Hz}$ and $7.1697\mathrm{Hz}$, respectively. The $\overline{f_s}$ of the UC in both cases is higher than the DC, which agrees with the discussion in Sumner et al. [14].

For the staggered arrangement considered in this study, the results for the general behaviour of the flow at $Re=3900$, in terms of flow patterns, are not dissimilar to the results in Sumner et al. [14], at $850<Re<1900$. The reason for this similarity is most likely due to both studies being in the lower range of the subcritical regime.

C18 is characterized by the highest values of $\overline{C_d}$ and the lowest values of $St$ among all nine flow patterns and where $\overline{C_{d_{DC}}}>\overline{C_{d_{UC}}}$. In C23, at the highest $\frac{P}{D}$ in the present study, the $St$ is close to the value for the single cylinder case. C25 and C26 are the only flow patterns in which the “flopping behaviour” can be observed with different values of $St$ for each cylinder. The difference between $S{t}_{UC}$ and $S{t}_{DC}$ decreased as the interference effects of close proximity decreased due to an increase in the distance between UC and DC.

4. Conclusions

Two-dimensional simulations of the flow around a single cylinder and two cylinders in various configurations were created using CFD software STAR-CCM+, version 12.06.011. The rectangular fluid domain was created with the dimensions $50D\times 25D$, where $D=0.01 m$. A mesh convergence analysis was conducted for the single cylinder at $Re=3900$ and the chosen mesh quality was then used for all simulations.

The initial simulations were for flow around a single cylinder at $Re=20, 40, 80, 200$ and $3900$. The drag coefficient and Strouhal numbers were obtained from these simulations and used to validate the computer model.

Two cylinders were then considered in tandem arrangement, at $\frac{L}{D}=1.5,3$ and $4$ for $Re=200$ and $Re=3900$. The results at $Re=200$ were found to be in good agreement with Meneghini et al. [6], providing further validation of the model. The results at $Re=3900$ showed that tandem arrangement is sensitive to $Re$ changes. The flow regime at $\frac{L}{D}=1.5$ was the “reattachment” regime, and transition to the “co-shedding” regime was observed at $\frac{L}{D}=3$ and $4$. Additionally, no lower gap limit was found for the “reattachment” indicating that the “extended body” regime does not exist at $Re=3900$. In the “reattachment” regime, at $Re=3900$, $St$ increases as the gap between cylinders decreases.

Two cylinders in the side-by-side arrangement were then presented, at $\frac{T}{D}=1.5,2$ and $3$ for $Re=200$ and $Re=3900$. The results in terms of flow regimes showed that the side-by-side arrangement is not as sensitive to $Re$ changes as the tandem arrangement. The “biased flow” regime at $\frac{T}{D}=1.5$ and $2$, and the “parallel vortex streets” regime at $\frac{T}{D}=3$, occurred for both $Re=200$ and $Re=3900$. Instabilities in the fluctuating forces observed in the “biased flow” provided evidence of flopping behaviour in the wake, which has been observed experimentally [29]. Vortex shedding frequencies were found to be different for the HC and the LC in the “biased flow”. The $\overline{f_s}$ values of both cylinders at $\frac{T}{D}=2$ are higher than at $\frac{T}{D}=1.5$, when $Re=200$. However, the $\overline{f_s}$ for both cylinders at $\frac{T}{D}=2$ are lower than at $\frac{T}{D}=1.5$ when $Re=3900$. This shows a level of sensitivity and differences between $Re$ and the gap ranges for each regime in the side-by-side configuration. In the results at $\frac{T}{D}=3$, the interference effects between the cylinders are not as strong as in the “biased flow”. The flopping behaviour was no longer present, and the magnitude of the forces on the cylinders was tending towards the single-cylinder case; however, the wakes of the HC and the LC remained synchronized in an anti-phase shedding pattern. The results in terms of flow patterns were in fairly good agreement with the work performed by Meneghini et al. [6], at $Re=200$. However, in terms of fluctuating forces, differences were observed. In particular, differences were seen in the values of $\overline{C_d}$ and $\overline{C_l}$ for $\frac{T}{D}=1.5$ and $\frac{T}{D}=2$ at $Re=200$, in contrast to the result of Meneghini et al. [6] where the same for HC and LC was reported in both cases. At $Re=3900$, larger differences between the averaged coefficients of the two cylinders indicate a dependence on $Re$as well as the gap.

The final simulations are related to the staggered arrangement at $Re=3900$. The flow patterns were grouped into three categories to compare them with the experimental study performed by Sumner et al. [14], in the range of $850<Re<1900$. The results showed that in terms of flow patterns and general behaviour of the flow, they are similar to the results at lower Reynolds numbers. This may be due to both studies being in the lower range of the subcritical regime. The key findings in the staggered configuration were that in the first group of flow patterns, at $\frac{P}{D}=1.0,\alpha =30°$, the highest value of $\overline{C_d}$ was found, along with the lowest value of $St$, and this was the only configuration where $\overline{C_{d_{DC}}}>\overline{C_{d_{UC}}}$. In the second group of flow patterns, at $\frac{P}{D}=4.0,\alpha =10°$, the cylinders are so far apart that the $St$ is close to the single cylinder case. In the third group of flow patterns, at $\frac{P}{D}=1.5,\alpha =60°$ and $\frac{P}{D}=2.5,\alpha =50°$, flopping behaviour was observed along with different values of $St$ for each cylinder. The $\overline{f_s}$ of the UC in both cases is higher than the DC.