# Numerical Investigation of Fluid Flow past Four Cylinders at Low Reynolds Numbers

Authors: Shafee Ahmad; Shams Ul-Islam
Journal: Mathematical Problems in Engineering
 (2021)

---

## Abstract

A numerical investigation on the effects of separation ratios and Reynolds numbers on the flow around four square cylinders in diamond arrangement has been carried out using the lattice Boltzmann method. The separation ratios between the cylinders vary fromg∗=1 to 15. The Reynolds numbers based on the diameter of the square cylinder and the inlet uniform inflow velocity are selected from Re = 80 to 160. The computations show that a total of five different flow regimes are observed over the selected ranges: single bluff-body, quasi-unsteady, chaotic flow, in-phase synchronized vortex shedding, and antiphase synchronized vortex shedding flow regimes. It is found that the flow features significantly depend on both the separation ratio and Reynolds number, with the former’s influence being more than the latter’s. We found that the critical spacing for four square cylinders in diamond arrangement for selected Reynolds numbers (80 ≤ Re ≤ 160) is in the range of 2 ≤ g∗ ≤ 5. The results reveal that the presence of secondary cylinder interaction frequencies indicates that, for chaotic flow regime, the wake pattern is not stable and there is a strong interaction of gap flows and continuous change in the direction of shed vortices behind the cylinders. The effects of the g∗ and Re on fluid forces, vortex shedding frequency, and flow separation have been examined in detail.


---

## Body

## 1. Introduction

Fluid-structure interactions occur in a number of engineering fields. Study on wake structure and fluid forces on multiple bluff bodies is of practical importance in a number of engineering applications, such as chimneys, offshore structures, towers, bridges, marine risers, and ocean platforms. It is also of importance in the analysis of fluid/structure interactions since the mechanism associated with the complex physical phenomena such as vortex shedding, flow separation, merging of vortices, and interaction of gap flows with the shed vortices is still not fully understood.A typical irregular variation of wake structures is when multiple bluff bodies is subject to a cross flow. When the separation ratiog∗ is too small between the cylinders, the vortices strongly interact, thus substantially increasing the drag forces. Separation ratio g∗ is defined as g∗=S/D, where S and D are the surface-to-surface distance between the cylinders and diameter of the cylinder, respectively. Investigations into the effects of the g∗ and Reynolds number (Re) on multiple bluff bodies have also been carried out. The Reynolds number is defined based on the diameter of the cylinder as Re=umaxD/v, where umax is the maximum velocity at the inlet and υ is the fluid kinematic viscosity. Kolar et al. [1] carried out an experimental investigation using the two-component laser Doppler velocimetry and discussed in detail the synchronized flow around two side-by-side square cylinders. The Reynolds number and separation ratio were focused on 23,100 and 2, respectively. Williamson [2] experimentally examined the flow features around two side-by-side circular cylinders at 0.5 < g∗ ≤ 5 at low Reynolds numbers varying from 50 to 200. According to Guillaume and LaRue [3], flow around two side-by-side cylinders experiences different drag values and vortex shedding frequencies. Kang [4] have carried out a detailed numerical examination of flow around two side-by-side circular cylinders. The author found different flow regimes while changing the Reynolds numbers from 40 to 160 with g∗<5. Agrawal et al. [5] conducted numerical investigation into the effects of g∗ on the flow past two side-by-side square cylinders at Re = 73 using the lattice Boltzmann method (LBM). A basic study of flow around two side-by-side square cylinders was carried out by Inoue et al. (2006). They studied the effects of the separation ratios and observed various flow regimes at Re = 150. Alam and Zhou [6] conducted an experimental study to examine the flow features around two side-by-side square cylinders, and they observed different flow regimes at Re = 300 and 1 ≤ g∗ ≤ 5. There has been much numerical work on the flow past multiple cylinders arranged side-by-side, such as the work by Ding et al. [7], Burattini and Agrawal [8], Han et al. [9], and Islam et al. [10]. It is known from these studies that the flow regimes strongly depend on the g∗ value between the cylinders.The separation ratio between the inline cylinders has great impact on the flow regimes and force statistics (Zdravkovich [11]). The experimental investigations carried out by Yen et al. [12] were for the Re and g∗ for flow around two tandem square cylinders. They observed different flow regimes by changing the values of the Re and g∗. Etminan et al. [13] did a finite volume study on flow past tandem square cylinders using finite volume method at low Reynolds numbers at g∗=5. The authors found that the flow changes its state from steady to unsteady at Re = 40. A number of numerical studies on flow past two tandem square cylinders have also been carried out with the motivation to examine the transition of flow regimes from one state to another by changing the spacing between the cylinders (see Sohankar [14] and Zhao et al. [15]).Experimental studies have been carried out for flow around four circular cylinders in an inline square configuration by Sayers [16]. Lam and Lo [17] experimentally examined various types of flow regimes and their corresponding Strouhal number (St) for four circular cylinders in an inline square configuration at Re = 2100 and g∗ ranging from 1.28 to 5.96. The authors observed different flow regimes by changing the value of g∗. Strouhal number St is defined as St=fsD/umax, where fs is the vortex shedding frequency. Farrant et al. [18] numerically examined vortex shedding from arrays of cylinders in various arrangements for different orientations and spacings at Re = 200 using the cell boundary element method. They found well-known flow features such as in-phase vortex shedding, antiphase vortex shedding, and vortex shedding flow regimes, which are well-known flow characteristics and observed mostly in case of two-cylinder arrangement. Song et al. [19] studied numerically various flow regimes and force statistics for flow around four square cylinders in an inline square configuration using the finite volume method (FVM). They systematically investigated the effect of spacing between the cylinders for Re = 300. The authors found that g∗=3.5 is the critical spacing and the flow regime changes from a bistable shielded flow regime to a vortex shedding flow regime by changing the values of g∗ from 1.5 to 8. Lam et al. [20] numerically studied flow interference effects of four cylinders at Re = 100 and 200 with 1.5 ≤ g∗ ≤ 5 and observed different flow regimes. Lin et al. [21] numerically studied the different flow regimes at Re = 200 for flow around four circular cylinders by changing the value of g∗. They found that the flow becomes critical at g∗=3. Abbasi et al. [22] numerically examined the effects of Re and g∗ on flow around four square cylinders in an inline square configuration. Experimental studies have been carried out for flow around four square cylinders in a square configuration (Liu et al. [23]). The authors experimentally measured the lift CL and drag CD forces at Re = 4.58 × 104 and g∗=3.45, 4.14, and 5.17 for four different angles of attacks. The authors found that the four cylinders forces are slightly different by changing the value of g∗. Lift coefficient is defined as CL=2Fy/ρumax2D, where Fy is the lift force in the transverse direction. Similarly, the drag coefficient is defined as CD=2Fx/ρumax2D, where Fx is the drag force in the streamwise direction. Zhang et al. [24] experimentally investigated the flow features and forces for array of four square cylinders at different g∗ and orientations at Re = 8000 using the Particle Image Velocimetry (PIV). They observed three different flow regimes while changing the values of pitch ratio from 2 to 5 and incidence angles from 0° to 45°. Gao et al. [25] performed numerical simulations of vortex-induced vibrations of four circular cylinders in a square configuration at a low Reynolds number of 150 and mass ratio of 2 using the commercial software package Ansys Fluent 16.0. They examined in detail the effect of incident angle α varying from 0° to 45° and reduced velocity Vr varying from 3 to 14. The authors found that the incident angle has a considerable effect on the vibration amplitudes and force coefficients of the four circular cylinders, especially for the cylinders arranged at the downstream positions, where the maximum force coefficients and vibration amplitudes occur. The fluid flow plays an important role in engineering applications [26]. The nature of the flow around multiple bluff bodies is highly dependent on the separation ratio between the cylinders and practically very important [27–30].Knowledge of the effect ofg∗ on the flow characteristics around the square cylinders in diamond arrangement is scarce. It is not clear how g∗ and Re affect the fluid dynamics as compared to isolated cylinder, tandemly arranged cylinders, and side-by-side arranged cylinders. Numerical and experimental studies for flow past four square cylinders in diamond arrangement are far from complete as compared to circular cylinders. More specifically, an example is how a change in g∗ from 1 to 15 will affect forces, flow transitions, and Strouhal numbers of the four cylinders, for a given Re. The main motivation of present numerical examination is to explore systematically further insight about the flow mechanism past four square cylinders in diamond arrangement, to get in-depth knowledge of important parameters, such as CD, CL, St, transition of flow regimes, and power spectra of CL by changing the values of Re and g∗. Furthermore, two important aspects will be presented in more detail: (i) merging of jet flows and shear layer reattachment at low spacings and (ii) variation in the wake size. Five Reynolds numbers and eleven separation ratios are investigated. Due to two-dimensional limitations, we chose Re up to 160 [31].The paper is organized as follows. In Section2, the details of lattice Boltzmann method, problem description, effect of computational domain, grid independence, and code validation study are presented. The computed results are presented and analyzed in Section 3. Finally, the most important conclusions are drawn in Section 4.


## 2. Numerical Method and Computational Details

Flow past four square cylinders in diamond shape is governed by the unsteady two-dimensional (2D) incompressible Navier–Stokes equations. The Cartesian coordinate system for the proposed problem is as follows:(1)∂ui∂t+uj∂ui∂xj=−1ρ∂p∂xi+v∂2ui∂xjxj,(2)∂ui∂xi=0.Here,p is the pressure, ρ is the fluid density, v is the kinematic viscosity of fluid, and ui (or uj) is the velocity component in the xi (or xj) coordinate direction.Instead of trying to discretize (1) and (2), a discretization of the Boltzmann equation is good enough to deal with the fluid flow problems. Computational fluid dynamics (CFD) has become very popular and plays an important role in analysis of fluid flow studies. The LBM has developed quickly in many aspects in the last twenty years and nowadays has become a powerful and attractive CFD tool for fluid flows. The LBM has some advantages compared to other CFD techniques. First, in the LBM the pressure is calculated from the equation of state [32]. This means once the density field is known, one can directly obtain the pressure fields. Second, no-slip boundary condition can be handled easily by bounce-back scheme [33]. Third, the fundamental Boltzmann equation is directly discretized [34]. Today, researchers use LBM in a number of disciplines [35–39]. In the above-mentioned references, one can easily find the streaming; collision; equilibrium distribution; discretization in time, space, and velocity; derivation of weighting coefficients; and recovery of the governing continuity and Navier–Stokes equations.The LBM can be derived from the Bhatnagar, Gross, and Krook (BGK) (1954) approximation of the Boltzmann equation [40] without body force:(3)∂tf+ζΔf=Ω.Here,fx,ζ,t is the particle distribution function in the phase space x,ζ. x and ζ are the position vector and the microscopic velocity, respectively. Ω is the collision operator. The simplest collision operator that can be used for Navier–Stokes simulations is the Bhatnagar–Gross–Krook (BGK) operator defined as follows:(4)Ω=feq−fτ.Here,feqx,ζ is the Maxwell–Boltzmann distribution function. τ is the single-relaxation time parameter which controls the stability. It is important to state here that the viscosity must be positive for the chosen LBM model which requires the choice of τ > 0.5. Now putting (4) into (3), we get(5)∂tf+ζ⋅Δf=feq−fτ.In the LBM-BGK method, for fluid flow, the discretized distribution functionfi is used. The discretized lattice Boltzmann equation [40] is(6)fix+eiδt,t+δt=fix,t+δtfieqx,t−fix,tτ.Here,fix,t is the density distribution function and is linked to the discrete velocity direction i (see Figure 1). ei, t, and δt are the discrete velocities, time, and time increment, respectively. The relaxation parameter τ is related to the fluid kinematic viscosity by v=cs2τ−0.5δt, where cs is the speed of sound. To determine the macroscopic equations for the LBM, one may need to perform a Taylor expansion of (6). Equation (6) can be divided into two main steps. The first part is collision (or relaxation):(7)fi∗x,t=fix,t+δtfieqx,t−fix,tτ.Figure 1
Lattice arrangement for the two-dimensional proposed problem.Here,fi∗ represents the distribution function after collisions. It is to be noted that collision is simply an algebraic local operation. The second part is propagation (or streaming):(8)fix+eiδt,t+δt=fi∗x,t.Thefieqx,t in (6) can be calculated as [41](9)fieqx,t=wiρ1+ei⋅ucs2+ei⋅u22cs4−u22cs2.Here,u is the macroscopic velocity vector and wis are the weighting coefficients. It is important to mention here that in lattice Boltzmann equation we can deal with the low Mach number Ma=δtumax/δxcs < 0.3 for near incompressible flows; therefore, the terms higher than u2 are ignored in the equilibria. d2q9 discrete velocity model is selected for this study (Figure 1). d and q represent the dimensions and number of particles, respectively. c2s is defined as c2/3 for d2q9 model. Here, c is the lattice speed and is defined as c=δx/δt. In this study, we choose lattice spacing δx as 1 lattice unit and time step δt as 1 ts. The magnitudes of weight coefficients wi are defined as(10)wi=49,i=0,19,i=1,2,3,4,136,i=5,6,7,8.The discrete velocities for thed2q9 model are given by [33](11)ei=0,i=0,ei=cosπ2i−1,sinπ2i−1,i=1,2,3,4,ei=2×cosπ2i−92,sinπ2i−92,i=5,6,7,8.Theρ and u can be obtained from(12)ρ=∑i=08fi=∑i=08fieq,u=1ρ∑i=08fi⋅ei=1ρ∑i=08fieq⋅ei.One can use the equation of state to calculate the pressure as given in the following:(13)p=ρcs2.From the Chapman–Enskog expansion [35, 39], one can easily recover the governing Navier–Stokes equations (1) and continuity (2) from the lattice Boltzmann method.In 2D case, the vorticity is calculated as(14)ωz=∂v∂x−∂u∂y.Figure2 shows a schematic view of four identical square cylinders in diamond arrangement, and the system is exposed in a coming uniform inflow velocity umax, where S is the spacing between the cylinders and D the diameter of the cylinder. The coming maximum velocity is assumed to be 2D. C1, C2, C3, and C4 represent the first, second, third, and fourth cylinders, respectively. The upstream distance (Lu) of the domain (the distance between the inlet boundary and the front surface of the first cylinder) and the downstream distance (Ld) of the computational domain (the distance between the rear surface of the second cylinder and the exit boundary) are assigned constant values of 10D and 35D, respectively. Lx and Ly are the streamwise and transverse directions, respectively. A Cartesian coordinate system (x, y) is used to describe the flow where x and y are incoming flow direction and perpendicular to the stream direction, respectively. CD and CL are the drag and lift forces in the streamwise and transverse directions, respectively. The variations of grid points for different separation values are presented in Table 1.Figure 2
Schematic diagram of the proposed problem.Table 1
Selected separation ratios for diamond-shaped arranged cylinders.
CasesLx × Lyg∗=1301 × 961g∗=2321 × 981g∗=3341 × 1001g∗=4361 × 1021g∗=5381 × 1041g∗=6401 × 1061g∗=7421 × 1081g∗=8441 × 1101g∗=9461 × 1121g∗=10481 × 1141g∗=15581 × 1241No-slip boundary conditions [33] are applied on the surfaces of the cylinders. Periodic boundary conditions [36] are applied at the top and bottom walls of the computational domain. Convective boundary condition [42] is incorporated at the exit of the domain. The forces on the surfaces of the cylinders are incorporated using the momentum exchange method [43]. The in-house computational code is developed for the Lahey Fortran Windows platform with 64 bits. All the computations are carried out in Intel i7-2600, dual core 3.2 GHz, 8 GB memory, and 320 GB hard disk for 64-bit Windows operating system. A typical calculation (10,00000-time steps) takes about 20 hours of CPU time. In Tables 2–4, the flow quantities mean drag coefficient (CDmean), St, root-mean-square values of drag coefficient (CDrms), and lift coefficient (CLrms) are presented. In Table 3, β = Ly/D is the blockage ratio of the computational domain. In Tables 2–5, the number within brackets gives the percentage deviation relative to the D = 20, β = 13, Lu = 10D, and Ld = 35D, respectively. In Tables 2–4, the parentage deviation in brackets is calculated from the following equation:(15)%error=absolutehighest value−lowest vauehighest value∗100.Table 2
Effect of cylinder size on flow past single square cylinder.
CDmeanStCDrmsCLrmsD = 101.5833 (3.9%)0.2068 (21%)0.0256 (20%)0.3593 (16%)D = 201.52210.16440.02060.3029D = 301.5108 (0.7%)0.1644 (0%)0.0198 (3.9%)0.3012 (0.6%)Table 3
Effect ofβ on flow past single square cylinder for Lu = 10D and Ld = 35D.
βCDmeanStCDrmsCLrms71.8208 (16.4%)0.1901 (13.5%)0.0237 (13%)0.3595 (15.7%)131.52210.16440.02060.3029191.4627 (4%)0.1601 (2.6%)0.0202 (1.9%)0.2978 (1.7%)Table 4
Influence ofLu on flow past isolated cylinder for β = 13 and Ld = 35D.
CasesCDmeanStCDrmsCLrmsLu = 5D1.5334 (0.7%)0.1671 (1.6%)0.0202 (1.9%)0.2911 (3.9%)Lu = 10D1.52210.16440.02060.3029Lu = 15D1.5381 (1%)0.1644 (0%)0.0211 (2.3%)0.3102 (2.2%)Table 5
Influence ofLd on flow past isolated cylinder for β = 13 and Lu = 10D.
CasesCDmeanStCDrmsCLrmsLd = 20D1.5229 (0.1%)0.1644 (0%)0.0204 (1%)0.3037 (0.3%)Ld = 35D1.52210.16440.02060.3029Ld = 50D1.5216 (0.03%)0.1644 (0%)0.0207 (0.5%)0.3030 (0.03%)Before starting simulations for proposed problem, the grid independence study is carried out first to select such grid size which does not affect the flow features and force statistics considerably. For this purpose, we choose 10-point grid, 20-point grid, and 30-point grid along the surface of the single square cylinder, at Re = 160. The flow parameters likeCDmean, St, CDrms, and CLrms obtained from these simulations are given in Table 2. It can be seen that the 10-point grid considerably affected the results in terms of percentage difference as compared to the 20-point grid. Furthermore, the 30-point grid takes more computational time as compared to the 20-point grid, and the difference in terms of percentage is almost negligible. Therefore, we have chosen D = 20 on each surface of the cylinder. Similar grid size was chosen by [44].We also studied the effect of computational domain for flow past isolated cylinder at Re = 160. The computational domain significantly affects the flow parameters and fastens or slows the vortex shedding process [45]. Table 3 indicates that β = 13 is a good choice for height of computational domain. There is some difference in terms of percentage observed between β = 13 and 19 for Lu = 10D and Ld = 35D. Table 4 shows that by fixing β = 13 and Ld = 35D, respectively, it can be seen that after Lu = 10D the influence of upstream location almost vanishes. Interestingly, at Lu = 15D the physical parameters show slight increasing behavior. Furthermore, by fixing β = 13 and Lu = 10D, respectively, there is negligible difference between the computed results at Ld = 35D and 50D. Even the difference with Ld = 20D is also negligible. We will use Ld = 35D as downstream distance from the rear surface of the C2 to the exit of the computational domain (see Table 5), so that there is reasonably large enough space for flow past multiple cylinders.Calculations were first carried out for a single square cylinder in a cross flow with Re = 80, 100, 120, 140, and 160 to obtain a baseline for comparison. The results forCDmean and St are compared with available experimental and numerical results in Figures 3(a)–3(d). The present CDmean results are very close to the numerical results of Malekzadeh and Sohankar [31]. The present CDrms results are almost identical to the numerical results of Sen et al. [46].Figure 3
Variation of (a)CDmean, (b) St, (c) CDrms, and (d) CLrms as a function of Re.
(a)(b)(c)(d)TheCDmean, St, CDrms, and CLrms all show a good agreement with the reported data. Figure 4(a) shows the alternate generation of negative and positive shed vortices, i.e., the well-known Karman vortex street. The streamline visualization in Figure 4(b) further confirms the regular generation of shed vortices behind the cylinder. Figures 4(c) and 4(d) show the CD and CL force histories. It can be observed that the CD settles to a periodic behavior with a mean value around 1.5221. The CL also settles to a regular sinusoidal function after the onset of wake instability leads to vortex shedding. The power spectra analysis of CL is presented in Figure 4(e). The Strouhal number is obtained using the fast Fourier transform (FFT). The frequency of the CL is found to be 0.1644.Figure 4
(a) Instantaneous vorticity contour visualization. (b) Streamline visualization. (c, d) Force time history. (e) Power spectra analysis of theCL, Re = 160.
(a)(b)(c)(d)(e)


## 3. Results and Discussion

The calculations are carried out for flow around four cylinders arranged in diamond shape. The separation ratio changes from 1 to 15 while the Re varies from 80 to 160. 1, 2, 3, and 4 in subscript of flow quantities represent the first, second, third, and fourth cylinders, respectively. For comparison, we also add the single-cylinder (SC) data in Figures5(a)–5(d) and 8(a)–8(d). We assign the names to flow regimes on the basis of wake interactions, shear layer reattachment, gap flows, variation of forces, and spectra behavior of lift coefficients. The vorticity contour and streamline visualization for flow around cylinders at different combinations of Re and g∗ are plotted in Figures 9(a)–9(f) and 10(a)–10(f) , respectively. These graphs help us understand the wake structure mechanism near downstream of the cylinders and its effect on the downstream region. At g∗=1 and Re = 100, there is no flow through the gaps. In this case, alternate vortex shedding is composed of negative and positive vortex forms behind the last cylinder (Figure 9(a)). Similar flow features were also observed for other chosen Reynolds numbers at g∗=1. Similar to that of a flow around an isolated cylinder, a single vortex street is observed behind the cylinders. This flow behavior is named as single bluff-body flow regime (SBBFR). The corresponding time variations of CD and CL for SBBFR are presented in Figures 11(a) and 11(b). This graph clearly indicates that there is no modulation in the signals. The drag coefficient for the C1 is relatively smaller than those for the other three cylinders. The power spectra further confirm the regular alternate generation of shed vortices behind the cylinder (see Figures 12(a)–12(d)). The power spectra also confirm that the primary vortex shedding dominates the flow and the small one extra harmonic in C2 does not affect the flow behind the cylinders. The same shedding values are observed for the side-by-side cylinders (St3 = St4 = 0.0820). Kang [4] numerically suggests that the single bluff-body flow regime occurs for flow around two side-by-side cylinders at g∗≤0.4. At g∗=2 and Re = 80, there is a shear layer reattachment from C1 onto C2. The flow through gaps is also not affected too much by the near downstream vortices. The flow pattern behind cylinders is almost steady and shows some unsteadiness near exit of the computational domain (Figure 9(b)). That is why this flow regime is named as quasi-unsteady flow regime (QUFR). At (Re, g∗) = (140, 2), the combination of shear layer and strong jet effect through the gaps considerably changes the flow characteristics behind the cylinders (Figure 9(c)). The irregular variation of vortex shedding is because the gap flows are considerably stronger. At g∗=2, the shed vortices behind the cylinders do not travel in alternate fashion as they move downstream and merge, and distortion can be seen in the flow field. The flow features of each cylinder are significantly influenced by upstream cylinder due to small spacing between the cylinders. The vortex shedding that occurs from two side-by-side cylinders (C3 and C4) does not have strong influence on the forces and wake structures. One can see noticeable modulations in CD and CL for cylinders C1 and C2 as compared to the other two cylinders (Figures 11(e) and 11(f)). This type of flow regime is named as chaotic flow regime (CFR). Lin et al. [21] had similar numerical observation for flow around four circular cylinders in an inline square configuration at Re = 200 and g∗=3. From Figure 9(d), we observe that the fully developed shed vortices from C1 reattach to C2 while at the downstream of C2 these shed vortices are not regular as they move downstream. Due to two gap flows, the shed vortices behind the cylinders C3 and C4 show some irregularity as they shed, and move downstream and interact with the shed vortices of C2. Again, the flow shows chaotic nature at far downstream of the domain. The flow regime is named as CFR again. At g∗=5 (Figure 11(g)), we observe that the shed vortices within the spacing have a strong interaction with the cylinder C2. Moreover, at the downstream position of C2, the shed vortices travel in strong fashion. As a result, one can see higher drag value for C1 compared to C2. In antiphase synchronized vortex shedding flow regime (ASVSFR), at Re = 100 and g∗=15, the vortices can be clearly seen within the gap between C1 and C2 as well as in the downstream region of C2 (Figure 9(e)), because the flow finds enough space between C1 and C2 to roll up in the form of shed vortices. Near exit of the domain, they are interacting with each other. The modulation in the CD and CL of C1 and C2 is due to the shed vortices within the gap of cylinders. Figure 9(e) confirms that, in case of g∗=15, the four single-cylinder wakes are completely synchronized. The well-known von Karman vortex street is clearly found behind each cylinder at this large separation ratio between the cylinders. The minor modulation in the CD2 and CL2 occurs due to minor difference between the vortex shedding frequencies of C2. According to the spectrum analyses, the primary vortex shedding frequency dominates the flow. At this separation ratio, the shedding frequencies of the cylinders C3 and C4 are close to the isolated cylinder case for all Reynolds numbers. Figure 9(f) confirms that, in case of g∗=9 and Re = 120, the four single-cylinder wakes are completely synchronized. In addition, the shedding behind cylinders C3 and C4 is in-phase. According to the spectrum analyses, the primary vortex shedding frequency dominates the flow. At this separation ratio, the shedding frequencies of the cylinders C3 and C4 are close to the isolated cylinder case for all Reynolds numbers. That is why this flow regime is named as in-phase synchronized vortex shedding flow regime (ISVSFR). Furthermore, in ASVSFR and ISVSFR, the cylinders C3 and C4 confirm the antiphase and in-phase synchronized vortex shedding behind these two cylinders. Previous experimental study of Kolar et al. [1] for flow around two side-by-side square cylinders and numerical study of Farrant et al. [18] for flow around four circular cylinders in an inline square configuration also suggest the in-phase synchronized vortex shedding for large separation ratio. From these graphs, it can be observed that, along with g∗, the Reynolds numbers also have significant influence on the change if flow features around four cylinders. Such observations are also consistent with the available experimental studies of Williamson [2] for Re = 100 to 200 and Zhou et al. [47] for Re = 150 to 450 for flow past two side-by-side cylinders.Figure 5
Variation ofCDmean as a function of g∗ for (a) Re = 80, (b) Re = 100, (c) Re = 120, (d) Re = 140, and (e) Re = 160.
(a)(b)(c)(d)(e)Figure 6
Variation ofSt as a function of g∗ for (a) Re = 80, (b) Re = 100, (c) Re = 120, (d) Re = 140, and (e) Re = 160.
(a)(b)(c)(d)(e)Figure 7
Variation ofCDrms as a function of g∗ for (a) Re = 80, (b) Re = 100, (c) Re = 120, (d) Re = 140, and (e) Re = 160.
(a)(b)(c)(d)(e)Figure 8
Variation ofCLrms as a function of g∗ for (a) Re = 80, (b) Re = 100, (c) Re = 120, (d) Re = 140, and (e) Re = 160.
(a)(b)(c)(d)(e)Figure 9
Vorticity contour representation at different combinations of Re andg∗. (a) g∗=1.0andRe=100. (b) g∗=2.0andRe=80. (c) g∗=2.0andRe=140. (d) g∗=5.0andRe=160. (e) g∗=15.0andRe=100. (f) g∗=9.0andRe=120.
(a)(b)(c)(d)(e)(f)Figure 10
Streamline visualization at different combinations of Re andg∗. (a) g∗=1.0andRe=100. (b) g∗=2.0andRe=80. (c) g∗=2.0andRe=140. (d) g∗=5.0andRe=160. (e) g∗=15.0andRe=100. (f) g∗=9.0andRe=120.
(a)(b)(c)(d)(e)(f)Figure 11
Time history analysis ofCD and CL at different values of Re and g∗. (a) g∗=1.0andRe=100. (b) g∗=1.0andRe=100. (c) g∗=2.0andRe=80. (d) g∗=2.0andRe=80. (e) g∗=2.0andRe=140. (f) g∗=2.0andRe=140. (g) g∗=5.0andRe=160. (h) g∗=5.0andRe=160. (i) g∗=15.0andRe=100. (j) g∗=15.0andRe=100. (k) g∗=9.0andRe=120. (l) g∗=9.0andRe=120.
(a)(b)(c)(d)(e)(f)(g)(h)(i)(j)(k)(l)Figure 12
Power spectra analysis ofCL at various values of g∗ and Re. (a) g∗=1.0andRe=100. (b) g∗=1.0andRe=100. (c) g∗=1.0andRe=100. (d) g∗=1.0andRe=100. (e) g∗=2.0andRe=80. (f) g∗=2.0andRe=80. (g) g∗=2.0andRe=80. (h) g∗=2.0andRe=80. (i) g∗=2.0andRe=140. (j) g∗=2.0andRe=140. (k) g∗=2.0andRe=140. (l) g∗=2.0andRe=140. (m) g∗=5.0andRe=160. (n) g∗=5.0andRe=160. (o) g∗=5.0andRe=160. (p) g∗=5.0andRe=160. (q) g∗=15.0andRe=100. (r) g∗=15.0andRe=100. (s) g∗=15.0andRe=100. (t) g∗=15.0andRe=100. (u) g∗=9.0andRe=120. (v) g∗=9.0andRe=120. (w) g∗=9.0andRe=120. (x) g∗=9.0andRe=120.
(a)(b)(c)(d)(e)(f)(g)(h)(i)(j)(k)(l)(m)(n)(o)(p)(q)(r)(s)(t)(u)(v)(w)(x)As presented in Figure10(b), there is no vortex shedding in the gap between C1 and C2, but small recirculation zones are clearly seen behind cylinder C1. The gap flows can generate large lift on cylinder C4. As a result, one can see only two small eddies behind the cylinders (Figure 10(a)). As g∗ is increased to 2, there is no recirculation within the gaps and some unsteadiness is seen at far downstream of the domain. These characteristics belong to QUFR and are observed only for Re = 80 at this specific separation ratio (Figure 10(b)). At (Re, g∗) = (140, 2), the interaction of shed vortices and interference becomes very complex and strong within the cylinders due to narrow jet flows through the gaps. The wake of C1 significantly shrinks and narrows, because the gap flows from each side rush into the central region. The shear layer reattachment from C1 onto C2 generates a very complex region behind C2, C3, and C4 in terms of small and large eddies (Figure 10(c)). When g∗ is increased to 5, the effect of gap flows is slightly reduced and the recirculation observed behind each cylinder further confirms that the wake deviates either downward or upward and does not affect the shed vortices from C1. This recirculation confirms the irregular generation of vortices observed at far downstream of the computational domain (Figure 10(d)). As a result of such instability, one can see some more peaks in the power spectra together with the primary vortex shedding frequency (see Figures 12(m)–12(p)). We believe that the CFR for the four square cylinders in diamond arrangement lies between 3 and 5 for all chosen Reynolds numbers. At g∗=9 and 15, the streamline visualization shows that the shear layer of the C1 rolls up into a fully developed vortex and impinges onto the cylinder C2. Due to reasonably large spacing in presence of C1 and C2, the vortex shedding of side-by-side cylinders C3 and C4 becomes antiphase for g∗=9 (Figure 10(f)) and in-phase for g∗=15 (Figure 10(e)). As a result of large separation ratio, the vortex shedding frequency of the cylinders is almost close to that of an isolated cylinder for all chosen Reynolds numbers.Figures11(a)–11(l) show the time variation of drag and lift forces for different observed flow regimes at different values of Re and g∗. At g∗=1 and Re = 100, the time analysis of CD and CL of four cylinders for single bluff-body flow regime shows periodic behavior without any modulations (Figures 11(a) and 11(b)). Furthermore, at g∗=2 and Re = 80, the drag coefficient of cylinder C1 shows steady behavior for QUFR (Figure 11(c)). The lowest drag observed for C2, CD3, and CD4 shows antiphase behavior. The existence of a secondary cylinder interaction frequencies and its contribution in CFR can be clearly seen from the irregular variation of drag and lift forces presented in Figures 11(e)–11(h) for different combinations of Re and g∗. The amplitude of forces continuously changes due to interaction of gap flows with the shed vortices behind the cylinders. The time variation of lift coefficients is becoming periodic for the ASVSFR and ISVSFR regimes except for the cylinder C2. The CL2 is strongly modulated and strongly interacts with other cylinders shed vortices. On the other hand, the drag coefficients show some minor modulations (Figures 11(i)–11(l)). In Figures 11(i)–11(l), the time signal analysis of lift coefficients confirms the existence of only primary vortex shedding frequency. Only CD2 and CL2 show some small modulation. This ensures that at reasonably large separation ratios the role of secondary cylinder interaction frequencies almost disappears for all four cylinders. The time variation of lift coefficients is sinusoidal for the three cylinders (C1, C3, C4), with the amplitude for the cylinder C2 being larger than that of the C1, C3, and C4 cylinders. Furthermore, at low separation ratios, the drag coefficient in the cylinders C3 and C4 is higher than that in the isolated cylinder. Zhou et al. [47] also reported that, at g∗=2 and Re = 0.35–1.04 × 104, the drag coefficients in the two side-by-side cylinders case were considerably higher than that in the isolated cylinder case.To further examine the transition of flow regimes, FFT analyses of the lift coefficients are performed, and some representative cases are presented in Figures12(a)–12(x). It can be seen in Figures 12(a)–12(x) that the influences of primary vortex shedding frequency and secondary cylinder interaction frequencies change once the flow regime changes from one state to another. Kumar et al. [48] proposed the secondary cylinder interaction frequency concept for the first time for flow around row of square cylinders. In SBBFR, the primary vortex shedding indicates that the vortex pattern behind the cylinders moves in alternate fashion. Different vortex frequencies (Figures 12(a) and 12(b)) are observed for cylinders C1 and C2, and the same shedding frequencies (Figures 12(c) and 12(d)) are observed for C3 and C4. One small peak is also found for C1 and C2. In QUFR, for (Re, g∗) = (80, 2), there is no evidence of other harmonics in the spectra and the primary vortex shedding frequency dominates the flow (Figures 12(e)–12(h)). Almost the same shedding frequencies are observed for all cylinders with negligible difference. This is due to the unsteadiness nature of the flow near exit of the computational domain. To examine the frequency characteristics of the chaotic flow regime (CFR), two representative cases ((Re, g∗) = (140, 2) and (160, 5)) are presented in Figures 12(i)–12(p). The presence of secondary cylinder interaction frequencies indicates that, for CFR, the wake pattern is not stable and there is a strong interaction of gap flows and continuous change in the direction of shed vortices behind the cylinders. In this study, we found that the chaotic nature increased with increasing Reynolds numbers at g∗=2. The highest peak in Figures 12(i)–12(l) (St1 = 0.0735 and St2 = 0.0731) for inline cylinders (C1 and C2) and (St3 = 0.0735 and St4 = 0.0731) for side-by-side cylinders (C3 and C4) is the primary vortex shedding frequency. The power spectrum analysis at (Re, g∗) = (160, 5) is shown in Figures 12(m)–12(p). It is found again that there exist secondary cylinder interaction frequencies even at reasonable spacing between the cylinders. A Strouhal number of 0.2060, 0.2075, 0.1765, and 0.1758 corresponds to the primary vortex shedding frequency for C1, C2, C3, and C4, respectively. In the chaotic flow regime, due to some spacing between the cylinders, the gap flow effect becomes smaller and only one to two harmonics are seen in the spectra. This confirms the interaction of shed vortices behind the cylinders at far downstream of the domain. The secondary cylinder interaction frequencies together with the primary vortex shedding frequency are also observed by Abbasi et al. [22] for flow past four square cylinders in an inline square configuration. The spectra analysis of lift coefficients for (Re, g∗) = (100, 15) and (120, 9) for all four cylinders is presented in Figures 12(q)–12(x) for ASVSFR and ISVSFR. It is observed that there exists only primary vortex shedding frequency for both flow regimes due to large spacing between the cylinders. Only one extra small peak is observed for C2 which is the way of C1. Interestingly, the power spectra of C2 still show one small peak together with the primary vortex shedding frequency at large separation ratio (g∗=15) shown in Figures 12(r) and 12(v).As shown in Figures13(a)–13(f), the streamwise velocity distribution at different cross-sections for different values of Re and g∗ shows different characteristics. It is clear from these figures that the streamwise velocity at different cross-sections is considerably affected as the flow transition occurs. For example, in case of CFR (Figures 13(c)–13(d)), the irregular variation of velocities even at far downstream location confirms the mutual interaction of generated shed vortices along with the gap flows. It is observed that the downstream cylinder C2 generally experiences large streamwise and transverse amplitudes.Figure 13
Streamwise velocity distribution at various cross-sections for different combinations ofg∗ and Re. (a) g∗=1.0andRe=100. (b) g∗=2.0andRe=80. (c) g∗=2.0andRe=140. (d) g∗=5.0andRe=160. (e) g∗=15.0andRe=100. (f) g∗=9.0andRe=120.
(a)(b)(c)(d)(e)(f)The results ofCDmean, St, CDrms, and CLrms of four cylinders as a function of g∗ for the different Re values are presented in Figures 5(a)–5(d), 6(a)–6(d), 7(a)–7(d), and 8(a)–8(d), respectively, together with the results of isolated cylinder for comparison. It can be seen from Figures 5(a)–5(d) for all chosen Reynolds numbers that the tandem cylinders (C1 and C2) are showing smaller mean drag forces than the side-by-side ones. As mean drag increases, wake interference decreases. In general, the CDmean values of C1 and C2 are considerably less than CDmean values of isolated cylinder at various separation ratios for all chosen Reynolds numbers. This confirms that the presence of cylinder in the wake of other cylinder reduces the forces. This finding agrees well with the results of Zdravkovich [11] for flow past tandem cylinders. Furthermore, we observed from the figures that CDmean3 and CDmean4 decrease with the increase in g∗, except at Re = 80 and g∗=2. This is due to the flow transition from single bluff-body flow regime to quasi-unsteady flow regime. However, the opposite trend is observed for the CDmean2; i.e., it decreases first to g∗=6, 5, 4, 4, and 3 for Re = 80, 100, 120, 140, and 160, respectively; then starts to increase; and again decreases by increasing the values of g∗. The maximum average drag force can be observed at g∗=1 for all Re values for CDmean2. As the spacing ratio increases, the CDmean3 and CDmean4 values approach SC value for all Reynolds numbers. The CDmean2 values are mostly affected up to g∗=5 due to critical flow regime observed at 2 ≤ g∗ ≤ 5. It is observed that the drag values are almost the same for cylinders C3 and C4 at g∗≥1 for all chosen Reynolds numbers. It is to be noted that, in case of two side-by-side cylinders, Guillaume and LaRue [3] experienced different drag values. In case of four cylinders in an inline square configuration, researchers observed negative drag values for the downstream cylinders (Sayers [16] and Abbasi et al. [22]). We did not find negative drag values for the downstream cylinders in diamond arrangement of four square cylinders.At small spacings between the cylinders due to the complex interferences among the cylinders, the variations of the wakes appear to be more complicated than that of a flow past a single square cylinder. Theg∗ is found to have a strong effect on the wake structures and forces of the cylinder C2, rather than the cylinders C3 and C4. It is clear from Figures 6(a)–6(e) that Strouhal number of all four cylinders shows abrupt changes up to g∗=5 for different Reynolds numbers. Then, it shows increasing behavior and almost approaches the SC value. At 1 ≤ g∗ ≤ 5 and Re = 80, 100, 120, 140, and 160, we found multiple peaks in the power spectra because of secondary cylinder interaction frequencies. Further, at g∗≥6, we observed that generally the Strouhal number of tandem cylinders (C1 and C2) as well as side-by-side cylinders (C3 and C4) is almost the same, which confirms that the primary vortex shedding frequencies are dominant. We further observed different vortex shedding frequencies for cylinders C1 and C2 for all selected Reynolds numbers at g∗=3. This is due to the strong interaction of gap flows with the shed vortices behind the cylinders. The results indicate that when the spacing between the cylinders is reasonably large, the St approaches single-cylinder value. We found that the forces are slightly different for each individual cylinder by changing the value of g∗ at fixed Re. Liu et al. [23] had similar experimental observations at Re = 4.58 × 104 and g∗=3.45, 4.14, and 5.17 for flow around four square cylinders. Furthermore, the results show that the vortex shedding frequency considerably depends on the separation ratio, especially at 1 ≤ g∗ ≤ 6. For the antiphase and in-phase synchronized vortex shedding flow regimes, the vortex shedding frequency of cylinders is almost constant at 9 ≤ g∗ ≤ 15, very close to that (St = 0.1516, 0.1596, 0.1644, 0.1644, and 0.1644 for Re = 80, 100, 120, 140, and 160, respectively) in the isolated cylinder case. Interestingly, the vortex shedding frequency drastically increases for cylinders C1 and C2 and decreases for cylinders C3 and C4 with increasing separation ratio for the chaotic flow regime.The effects of separation ratiosg∗ on CDrms of four cylinders are shown in Figures 7(a)–7(e). It can be seen from Figures 7(a)–7(e) that with g∗ the values of CDrms are generally reduced and approach SC value except at CDrms2 for all chosen Reynolds numbers. At g∗=1, the CDrms1 values reach their minimum, and the other three cylinders attain the maximum values for all selected ranges of Reynolds numbers. It is observed that, due to chaotic flow regime, the effect of g∗ for all Reynolds numbers on the four cylinders CDrms values is significant and becomes small when g∗≥7. As g∗ increases, as shown in Figures 7(a)–7(e), the variation of CDrms1 is almost constant with minor variation in that of the SC. The CDrms1, CDrms3, and CDrms4 values appear to be smaller than those of the CDrms2, except for the g∗=1 and 2 for all chosen Reynolds numbers. At large separation ratio g∗=7 to 9, the CDrms of C3 and C4 is maintained at constant values, and the values for the cylinder C2 is higher than those for the cylinders C3 and C4. The differences of CDrms2 from CDrms3 and CDrms4 for g∗=7 to 9 become higher as the Reynolds number increased. When g∗ increases to 15, the CDrms1, CDrms3, and CDrms4 are almost the same. However, the variation of CDrms2 is quite different.In Figures8(a)–8(e), the calculation results of CLrms for various separation ratios at 80 ≤ Re ≤ 160 are compared with the isolated cylinder data. It can be seen that, due to the shear layer reattachment from C1 onto C2 and the complex interaction of gap flows with the shed vortices behind C3 and C4 from g∗=2 to 5, the variation of the CLrms with separation ratios for four cylinders appears to be more complicated than that of the SC value. The CLrms3 and CLrms4 decrease with increasing g∗ but are still larger than those of CLrms1. The results show that the CLrms2 is considerably affected by changing the value of g∗ at fixed Reynolds number. A possible explanation for such observation is that the wake structure behavior for the cylinder C2 is too complicated than that for the other cylinders, involving the gap flows from the upstream cylinders and the shear layers interaction with cylinder C1. It is seen from the CLrms values that the flow regimes are very sensitive to the separation ratio up to g∗≥7. It is found that CLrms1 value is lower than the isolated cylinder value for all Reynolds numbers up to g∗=7 and then approach SC value for g∗≥8. This mean that the effect of shear layer becomes less beyond g∗=7. The root-mean-square values of C3 and C4 cylinders are almost the same for all selected ranges and approach SC value. The drastic changes in CLrms2 by changing the values of g∗ for fixed Reynolds number confirm that the secondary cylinder interaction frequencies occur continuously due to gap flows interaction with shed vortices behind the cylinders. These characteristics further tell us that the shed vortices near downstream of C2 are affected by the vortices generated from C1. Furthermore, with increasing g∗, the CLrms1 value first decreases to g∗=3, then slightly increases, and finally approaches the SC value.The graphical representation of various flow regimes, found in this numerical study, is given in Figure14 at all separation ratios and Reynolds numbers. It can be seen that, at low separation ratios (g∗=1) and Re = 80, 100, 120, 140, and 150, the SBBFR occurs. In the small separation ratio (g∗=1), the gap flow has almost no effect, and the shear layer on the C1 cannot fall off freely. The four cylinders can act as a single bluff body, as shown in Figure 9(a). At (Re, g∗) = (80, 2), (100, 2), and (120, 2), the QUFR occurs. On the other hand, at (Re, g∗) = (80, 3–5), (100, 3–5), (120, 3–5), (140, 2–5), and (160, 2–5), the CFR is observed. At (Re, g∗) = (80, 6–10 and 15), (100, 6–10), (120, 6–9), (140, 7–10 and 15), and (160, 6 and 8–10), the ISVSFR can be clearly seen. For (Re, g∗) = (100, 15), (120, 10 and 15), (140, 6), and (160, 7 and 15), the ASVSFR occurs. The flow features vary only slightly as Re increases for g∗=6 to 15. The ISVSFR is dominant as compared to ASVSFR. As shown in Figure 14, the critical spacing for four square cylinders in diamond arrangement for selected Reynolds numbers (80 ≤ Re ≤ 160) is in the range of 2 ≤ g∗ ≤ 5. It is to be noted that for the four cylinders in diamond arrangement, the critical spacing appears to be generally wider than that of two side-by-side square cylinders (Agrawal et al. [5]) and two tandem cylinders (Yen et al. [12]). It can be seen that the in-phase and antiphase nature of the cylinders C3 and C4 occurred when the separation ratio increased and the fully developed shed vortices existed between the cylinders C1 and C2 at higher separation ratios. Williamson [2] also experimentally observed that the wakes behind the two side-by-side cylinders were either antiphase or in-phase when the separation ratio increased.Figure 14
Flow regimes as a function of Re andg∗.


## 4. Conclusions

Two-dimensional computations have been performed to examine the influence of Re andg∗ of four square cylinders in a diamond arrangement. The effects of the Re and g∗ on the vortex shedding modes and dynamic responses of the four square cylinders are examined. The main conclusions are summarized as follows:(i)
Numerical results showed that, by varyingg∗ from 1 to 15 and Re from 80 to 160, a total of five different kinds of flow regimes were identified: single bluff-body flow regime (80 ≤ Re ≤ 160 and g∗=1), quasi-unsteady flow regime (Re = 80, 100, and 120 and g∗=2), critical flow regime (80 ≤ Re ≤ 160 and 3 ≤ g∗ ≤ 5; also at g∗=2 for Re = 140 and 160), antiphase synchronized vortex shedding flow regime ((Re, g∗) = (140 and 160, 2), (80, 100, 120, and 160, 7), (80 ≤ Re ≤ 160, 8), (80 ≤ Re ≤ 160, 9), (80, 100, 140, and 160, 10), (80 and 140, 15)), and in-phase synchronized vortex shedding flow regime ((Re, g∗) = (140, 6), (160, 7), (120, 10), (100, 120, and 160, 15)). It was found that the flow regimes strongly depended on the separation ratio as well as the Reynolds number.(ii)
The vortex shedding frequency strongly depended on the separation ratio. The secondary cylinder interaction frequencies also play a strong role for the chaotic flow regime. It is also found that the modulations inCD and CL for cylinders C1 and C2 in chaotic flow regime are significant. With increasing separation ratio, the primary vortex shedding frequency dominates the flow and remains nearly constant for the antiphase and in-phase synchronized vortex shedding flow regimes. In these flow regimes, the flow finds enough space between C1 and C2 to roll up in the form of shed vortices. We found the well-known von Karman vortex street behind each cylinder at larger separation ratio between the cylinders. In case of quasi-unsteady flow regime, we found the steady behavior of drag force of cylinder C2.(iii)
The secondary cylinder interaction frequencies dominate the flow in case of critical flow regime. On the other hand, the primary vortex shedding frequency dominates the flow at larger separation ratios.(iv)
It is found that theCDmean3 and CDmean4 decrease with the increase in g∗, except at Re = 80 and g∗=2. This is due to the flow transition from single bluff-body flow regime to quasi-unsteady flow regime. However, the opposite trend is observed for the CDmean2; i.e., it decreases first to g∗=6, 5, 4, 4, and 3 for Re = 80, 100, 120, 140, and 160, respectively; then starts to increase; and again decreases by increasing the values of g∗. The maximum average drag force can be observed at g∗=1 for all Re values for CDmean2. As the spacing ratio increases, the CDmean3 and CDmean4 values approach the SC value for all Reynolds numbers. The CDmean2 value are mostly affected up to g∗=5 due to critical flow regime observed at 2 ≤ g∗ ≤ 5. We found that the vortex shedding frequencies are slightly different for each individual cylinder by changing the value of g∗ at fixed Re. Furthermore, we also found that the vortex shedding frequency considerably depends on the separation ratio, especially at 1 ≤ g∗ ≤ 6. For the antiphase and in-phase synchronized vortex shedding flow regimes, the vortex shedding frequency of cylinders is almost constant at 9 ≤ g∗ ≤ 15, very close to that (St = 0.1516, 0.1596, 0.1644, 0.1644, and 0.1644 for Re = 80, 100, 120, 140, and 160, respectively) in the isolated cylinder case. In general, the CDrms1 and CLrms1 have larger values compared to the other cylinders, C2, C3, and C4, for all chosen Re and g∗.



---

*Source: 1127324-2021-07-27.xml*