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

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. +e separation ratios between the cylinders vary from g∗ � 1 to 15. +e Reynolds numbers based on the diameter of the square cylinder and the inlet uniform inflow velocity are selected from Re� 80 to 160. +e 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. +e 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. +e effects of the g∗ and Re on fluid forces, vortex shedding frequency, and flow separation have been examined in detail.


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 ratio (g * ) 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. e Reynolds number is defined based on the diameter of the cylinder as Re � u max D/v, where u max 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 sideby-side square cylinders. e 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. e 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). ey 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. ere 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. e separation ratio between the inline cylinders has great impact on the flow regimes and force statistics (Zdravkovich [11]). e experimental investigations carried out by Yen et al. [12] were for the Re and g * for flow around two tandem square cylinders. ey 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. e 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. e authors observed different flow regimes by changing the value of g * . Strouhal number St is defined as St � f s D/u max , where f s 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. ey found well-known flow features such as in-phase vortex shedding, antiphase vortex shedding, and vortex shedding flow regimes, which are wellknown 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). ey systematically investigated the effect of spacing between the cylinders for Re � 300. e 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 * . ey 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]). e authors experimentally measured the lift (C L ) and drag (C D ) forces at Re � 4.58 × 10 4 and g * � 3.45, 4.14, and 5.17 for four different angles of attacks. e authors found that the four cylinders forces are slightly different by changing the value of g * . Lift coefficient is defined as where F y is the lift force in the transverse direction. Similarly, the drag coefficient is defined as where F x 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). ey 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 vortexinduced 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. ey examined in detail the effect of incident angle α varying from 0°to 45°and reduced velocity V r varying from 3 to 14. e 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. e fluid flow plays an important role in engineering applications [26]. e nature of the flow around multiple bluff bodies is highly dependent on the separation ratio between the cylinders and practically very important [27][28][29][30].
Knowledge of the effect of g * 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. e 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 C D , C L , St, transition of flow regimes, and power spectra of C L 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]. e paper is organized as follows. In Section 2, the details of lattice Boltzmann method, problem description, effect of computational domain, grid independence, and code validation study are presented.
e computed results are presented and analyzed in Section 3. Finally, the most important conclusions are drawn in Section 4.

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. e Cartesian coordinate system for the proposed problem is as follows: Here, p is the pressure, ρ is the fluid density, v is the kinematic viscosity of fluid, and u i (or u j ) is the velocity component in the x i (or x j ) 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. e 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. e LBM has some advantages compared to other CFD techniques. First, in the LBM the pressure is calculated from the equation of state [32]. is 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].
ird, the fundamental Boltzmann equation is directly discretized [34]. Today, researchers use LBM in a number of disciplines [35][36][37][38][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. e LBM can be derived from the Bhatnagar, Gross, and Krook (BGK) (1954) approximation of the Boltzmann equation [40] without body force: Here, f(x, ζ, 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. e simplest collision operator that can be used for Navier-Stokes simulations is the Bhatnagar-Gross-Krook (BGK) operator defined as follows: Here, f eq (x, ζ) 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 In the LBM-BGK method, for fluid flow, the discretized distribution function f i is used. e discretized lattice Boltzmann equation [40] is Here, f i (x, t) is the density distribution function and is linked to the discrete velocity direction i (see Figure 1). e i , t, and δt are the discrete velocities, time, and time increment, respectively. e relaxation parameter τ is related to the fluid kinematic viscosity by v � c 2 s (τ − 0.5)δt, where c s 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. e first part is collision (or relaxation): Here, f * i represents the distribution function after collisions. It is to be noted that collision is simply an algebraic local operation. e second part is propagation (or streaming): (6) can be calculated as [41] f Here, u is the macroscopic velocity vector and w i s are the weighting coefficients. It is important to mention here that in lattice Boltzmann equation we can deal with the low Mach number (Ma � δtu max /δxc s ) < 0.3 for near incompressible flows; therefore, the terms higher than u 2 are ignored in the equilibria. d 2 q 9 discrete velocity model is selected for this study (Figure 1). d and q represent the dimensions and number of particles, respectively. c 2 s is defined as c 2 /3 for d 2 q 9 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. e magnitudes of weight coefficients w i are defined as , i � 0, e discrete velocities for the d 2 q 9 model are given by [33] Mathematical Problems in Engineering e ρ and u can be obtained from One can use the equation of state to calculate the pressure as given in the following: 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 Figure 2 shows a schematic view of four identical square cylinders in diamond arrangement, and the system is exposed in a coming uniform inflow velocity u max , where S is the spacing between the cylinders and D the diameter of the cylinder. e coming maximum velocity is assumed to be 2D. C 1 , C 2 , C 3 , and C 4 represent the first, second, third, and fourth cylinders, respectively. e upstream distance (L u ) of the domain (the distance between the inlet boundary and the front surface of the first cylinder) and the downstream distance (L d ) 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. L x and L y 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. C D and C L are the drag and lift forces in the streamwise and transverse directions, respectively. e variations of grid points for different separation values are presented in Table 1.
No-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. e forces on the surfaces of the cylinders are incorporated using the momentum exchange method [43]. e 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 (C Dmean ), St, root-mean-square values of drag coefficient (C Drms ), and lift coefficient (C Lrms ) are presented. In Table 3, β � L y /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, L u � 10D, and L d � 35D, respectively. In Tables 2-4, the parentage deviation in brackets is calculated from the following equation: 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, 20point grid, and 30-point grid along the surface of the single square cylinder, at Re � 160. e flow parameters like C Dmean , St, C Drms , and C Lrms 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 20point grid, and the difference in terms of percentage is almost negligible. erefore, 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. e 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.
ere is some difference in terms of percentage observed between β � 13 and 19 for L u � 10D and L d � 35D. Table 4 shows that by fixing β � 13 and L d � 35D, respectively, it can be seen that after L u � 10D the influence of upstream location almost vanishes. Interestingly, at L u � 15D the physical parameters show slight increasing behavior. Furthermore, by fixing β � 13 and L u � 10D, respectively, there is negligible difference between the computed results at L d � 35D and 50D. Even the difference with L d � 20D is also negligible. We will use L d � 35D as downstream distance from the rear surface of the C 2 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. e results for C Dmean and St are compared with available experimental and numerical results in Figures 3(a)-3(d). e present C Dmean results are very close to the numerical results of Malekzadeh and Sohankar [31]. e present C Drms results are almost identical to the numerical results of Sen et al. [46]. Table 1: Selected separation ratios for diamond-shaped arranged cylinders.
Cases   Table 4: Influence of L u on flow past isolated cylinder for β � 13 and L d � 35D.
Cases Figure  and 4(d) show the C D and C L force histories. It can be observed that the C D settles to a periodic behavior with a mean value around 1.5221. e C L also settles to a regular sinusoidal function after the onset of wake instability leads to vortex shedding. e power spectra analysis of C L is presented in Figure 4(e). e Strouhal number is obtained using the fast Fourier transform (FFT). e frequency of the C L is found to be 0.1644.

Results and Discussion
e calculations are carried out for flow around four cylinders arranged in diamond shape. e 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 Figures 5(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. e 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. ese 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. is flow behavior is named as single bluff-body flow regime (SBBFR). e corresponding time variations of C D and C L for SBBFR are presented in Figures 11(a) and 11(b). is graph clearly indicates that there is no modulation in the signals. e drag coefficient for the C 1 is relatively smaller than those for the other three cylinders. e power spectra further confirm the regular alternate generation of shed vortices behind the cylinder (see Figures 12(a)-12(d)). e power spectra also confirm that the primary vortex shedding dominates the flow and the small one extra harmonic in C 2 does not affect the flow behind the cylinders. e 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 C 1 onto C 2 . e flow through gaps is also not affected too much by the near downstream vortices. e flow pattern behind cylinders is almost steady and shows some unsteadiness near exit of the computational domain (Figure 9(b)). at 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)). e 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.
e flow features of each cylinder are significantly influenced by upstream cylinder due to small spacing between the cylinders. e vortex shedding that occurs from two sideby-side cylinders (C 3 and C 4 ) does not have strong influence on the forces and wake structures. One can see noticeable modulations in C D and C L for cylinders C 1 and C 2 as compared to the other two cylinders (Figures 11(e) and 11(f)).
is 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 C 1 reattach to C 2 while at the downstream of C 2 these shed vortices are not regular as they move downstream. Due to two gap flows, the shed vortices behind the cylinders C 3 and C 4 show some irregularity as they shed, and move downstream and interact with the shed vortices of C 2 . Again, the flow shows chaotic nature at far downstream of the domain. e 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 C 2 . Moreover, at the downstream position of C 2 , the shed vortices travel in strong fashion. As a result, one can see higher drag value for C 1 compared to C 2 . In antiphase synchronized vortex shedding flow regime (ASVSFR), at Re � 100 and g * � 15, the vortices can be clearly seen within the gap between C 1 and C 2 as well as in the downstream region of C 2 (Figure 9(e)), because the flow finds enough space between C 1 and C 2 to roll up in the form of shed vortices. Near exit of the domain, they are interacting with each other. e modulation in the C D and C L of C 1 and C 2 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. e well-known von Karman vortex street is clearly found behind each cylinder at this large separation ratio between the cylinders. e minor modulation in the C D2 and C L2 occurs due to minor difference between the vortex shedding frequencies of C 2 . According to the spectrum analyses, the primary vortex shedding frequency dominates the flow. At this separation ratio, the shedding frequencies of the cylinders C 3 and C 4 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 C 3 and C 4 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 C 3 and C 4 are close to the isolated cylinder case for all Reynolds numbers. at is why this flow regime is named as in-phase synchronized vortex shedding flow regime (ISVSFR). Furthermore, in ASVSFR and ISVSFR, the cylinders C 3 and C 4 confirm the antiphase and in-  Mathematical Problems in Engineering 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.
As presented in Figure 10(b), there is no vortex shedding in the gap between C 1 and C 2 , but small recirculation zones are clearly seen behind cylinder C 1 . e gap flows can generate large lift on cylinder C 4 . 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.
ese 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. e wake of C 1 significantly shrinks and narrows, because the gap  Mathematical Problems in Engineering flows from each side rush into the central region. e shear layer reattachment from C 1 onto C 2 generates a very complex region behind C 2 , C 3 , and C 4 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 C 1 . is 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 C 1 rolls up into a fully developed vortex and impinges onto the cylinder C 2 . Due to reasonably large spacing in presence of C 1 and C 2 , the vortex shedding of side-by-side cylinders C 3 and C 4 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. Figures 11(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 C D and C L 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 C 1 shows steady behavior for QUFR (Figure 11(c)). e lowest drag observed for C 2 , C D3 , and C D4 shows antiphase behavior. e 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)- Only C D2 and C L2 show some small modulation. is ensures that at reasonably large separation ratios the role of secondary cylinder interaction frequencies almost disappears for all four cylinders. e time variation of lift coefficients is sinusoidal for the three cylinders (C 1 , C 3 , C 4 ), with the amplitude for the cylinder C 2 being larger than that of the C 1 , C 3 , and C 4 cylinders. Furthermore, at low separation ratios, the drag coefficient in the cylinders C 3 and C 4 is higher than that in the isolated cylinder. Zhou et al. [47] also reported that, at g * � 2 and Re � 0.35-1.04 × 10 4 , 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 Figures 12(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 C 1 and C 2 , and the same shedding frequencies (Figures 12(c) and 12(d)) are observed for C 3 and C 4 . One small peak is also found for C 1 and C 2 . 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. is 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). e 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. e highest peak in Figures 12(i)-12(l) (St1 � 0.0735 and St2 � 0.0731) for inline cylinders (C 1 and C 2 ) and (St3 � 0.0735 and St4 � 0.0731) for side-by-side cylinders (C 3 and C 4 ) is the primary vortex shedding frequency. e 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 C 1 , C 2 , C 3 , and C 4 , 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. is confirms the interaction of shed vortices behind the cylinders at far downstream of the domain. e 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. e spectra   Lift coefficients    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 C 2 which is the way of C 1 . Interestingly, the power spectra of C 2 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 Figures 13(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 crosssections 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 C 2 generally experiences large streamwise and transverse amplitudes.
e results of C Dmean , St, C Drms , and C Lrms 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 (C 1 and C 2 ) are showing smaller mean drag forces than the side-by-side ones. As mean drag increases, wake interference decreases. In general, the C Dmean values of C 1 and C 2 are considerably less than C Dmean values of isolated cylinder at various separation ratios for all chosen Reynolds numbers. is confirms that the presence of cylinder in the wake of other cylinder reduces the forces. is finding agrees well with the results of Zdravkovich [11] for flow past tandem cylinders. Furthermore, we observed from the figures that C Dmean3 and C Dmean4 decrease with the increase in g * , except at Re � 80 and g * � 2. is 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 C Dmean2 ; 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 * . e maximum average drag force can be observed at g * � 1 for all Re values for C Dmean2 . As the spacing ratio increases, the C Dmean3 and C Dmean4 values approach SC value for all Reynolds numbers. e C Dmean2 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 C 3 and C 4 at g * ≥ 1 for all chosen Reynolds numbers. It is to be noted that, in case of two side-byside 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. e g * is found to have a strong effect on the wake structures and forces of the cylinder C 2 , rather than the cylinders C 3 and C 4 . 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. en, 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 (C 1 and C 2 ) as well as side-by-side cylinders (C 3 and C 4 ) is almost the same, which confirms that the primary vortex shedding frequencies are dominant. We further observed different vortex shedding frequencies for cylinders C 1 and C 2 for all selected Reynolds numbers at g * � 3.
is is due to the strong interaction of gap flows with the shed vortices behind the cylinders.
e 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 × 10 4 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 inphase 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 C 1 and C 2 and decreases for cylinders C 3 and C 4 with increasing separation ratio for the chaotic flow regime. e effects of separation ratios g * on C Drms 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 C Drms are generally reduced and approach SC value except at C Drms2 for all chosen Reynolds numbers. At g * � 1, the C Drms1 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 C Drms values is significant and becomes small when g * ≥ 7. As g * increases, as shown in Figures 7(a)-7(e), the variation of C Drms1 is almost constant with minor variation in that of the SC. e C Drms1 , C Drms3 , and C Drms4 values appear to be smaller than those of the C Drms2 , except for the g * � 1 and 2 for all chosen Reynolds numbers. At large separation ratio g * � 7 to 9, the C Drms of C 3 and C 4 is maintained at constant values, and the values for the cylinder C 2 is higher than those for the cylinders C 3 and C 4 .
e differences of C Drms2 from C Drms3 and C Drms4 for g * � 7 to 9 become higher as the Reynolds number increased. When g * increases to 15, the C Drms1 , C Drms3 , and C Drms4 are almost the same. However, the variation of C Drms2 is quite different.
In Figures 8(a)-8(e), the calculation results of C Lrms 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 C 1 onto C 2 and the complex interaction of gap flows with the shed vortices behind C 3 and C 4 from g * � 2 to 5, the variation of the C Lrms with separation ratios for four cylinders appears to be more complicated than that of the SC value. e C Lrms3 and C Lrms4 decrease with increasing g * but are still larger than those of C Lrms1 . e results show that the C Lrms2 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 C 2 is too complicated than that for the other cylinders, involving the gap flows from the upstream cylinders and the shear layers interaction with cylinder C 1 . It is seen from the C Lrms values that the flow regimes are very sensitive to the separation ratio up to g * ≥ 7. It is found that C Lrms1 value is lower than the isolated cylinder value for all Reynolds numbers up to g * � 7 and then approach SC value for g * ≥ 8. is mean that the effect of shear layer becomes less beyond g * � 7. e root-mean-square values of C 3 and C 4 cylinders are almost the same for all selected ranges and approach SC value. e drastic changes in C Lrms2 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. ese characteristics further tell us that the shed vortices near downstream of C 2 are affected by the vortices generated from C 1 . Furthermore, with increasing g * , the C Lrms1 value first decreases to g * � 3, then slightly increases, and finally approaches the SC value. e graphical representation of various flow regimes, found in this numerical study, is given in Figure 14 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 C 1 cannot fall off freely. e 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. e flow features vary only slightly as Re increases for g * � 6 to 15. e 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-byside 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 C 3 and C 4 occurred when the separation ratio increased and the fully developed shed vortices existed between the cylinders C 1 and C 2 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.
(ii) e vortex shedding frequency strongly depended on the separation ratio. e secondary cylinder interaction frequencies also play a strong role for the chaotic flow regime. It is also found that the modulations in C D and C L for cylinders C 1 and C 2 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 C 1 and C 2 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 C 2 .
(iii) e 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 the C Dmean3 and C Dmean4 decrease with the increase in g * , except at Re � 80 and g * � 2. is 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 C Dmean2 ; 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 * . e maximum average drag force can be observed at g * � 1 for all Re values for C Dmean2 . As the spacing ratio increases, the C Dmean3 and C Dmean4 values approach the SC value for all Reynolds numbers.
e C Dmean2 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 C Drms1 and C Lrms1 have larger values compared to the other cylinders, C 2 , C 3 , and C 4 , for all chosen Re and g * .

Nomenclature
C D � 2F x /ρu 2 max D: Drag coefficient C L � 2F y /ρu 2 max D: Lift coefficient C Dmean : Mean drag coefficient C Drms : Root-mean-square value of drag coefficient C Lrms : Root-mean-square value of lift coefficient D: Size of the cylinder (m) e i : Discrete velocities f s : Vortex shedding frequency f: Particle distribution function f (eq) : Maxwell-Boltzmann distribution function F x : Drag force in the streamwise direction F y : Lift force in the transverse direction g * : Separation ratio L u : Upstream distance from the inlet position L d : Downstream distance up to the exit of the computational domain L x : Computational length L y : Computational height Ma: Mach number p: Pressure Re � u max D/υ:

Reynolds number S:
Surface-to-surface distance between the cylinders St � f s D/u max : Strouhal number u max : Maximum inflow velocity u: Flow speed (m/s) w i : Weighting coefficients β: Blockage ratio ρ: Density of the fluid (kg/m 3 ) υ: Kinematic viscosity of the fluid (m 2 /s) τ: Single-relaxation-time parameter Ω: Collision operator.
Data Availability e data that support the main findings of this numerical study are available from the corresponding author upon request.

Conflicts of Interest
e authors declare that there are no conflicts of interest regarding the submission and publication of this manuscript.