Transitions in the flow patterns and aerodynamic characteristics of the flow around staggered rows of cylinders

A two-dimensional numerical study of flow across rows of identical square cylinders arranged in staggered fashion is carried out. This study will unreveal complex flow physics depending upon the Reynolds number (Re) and gap spacing (g) between the cylinders. The combined effect of Reynolds number and gap spacing on the flow physics around staggered rows of cylinders are numerically studied for 20 ≤ Re ≤ 140 and 1 ≤ g ≤ 6. We use the lattice Boltzmann method for numerical computations. It is found that with increase in gap spacing between the cylinders the critical Reynolds number for the onset of vortex shedding also increases. We observed a strong effect of Reynolds number at g = 2 and 4. Secondary cylinder interaction frequency disappears for large Reynolds number at g = 6 and 5 and the flow around cylinders are fully dominated by the primary vortex shedding frequency. This ensures that at large gap spacing with an increase in the Reynolds number the wakes interaction between and behind the cylinders is weaken. Furthermore, it also ensures that the wake interaction behind the cylinders is strongly influenced by the jets in the gap spacing between the cylinders. We also found that g = 2 is the critical gap spacing for flow across rows of staggered square cylinders for the considered range of Reynolds number. Depending on the Reynolds number we observed; synchronous, quasi-periodic-I, quasi-periodic-II, and chaotic flow patterns. In synchronous flow pattern, an in-phase and anti-phase characteristics of consecutive cylinders has been observed. The important physical parameters are also analyzed and discussed in detail.


Introduction
The complex flow physics in terms of wake structure mechanism and formation of two rows of staggered identical square cylinders has attracted researchers' attention due to its importance in understanding industrial flows, for example, overhead power-line bundles, tubes in heat exchangers, chimney stacks, stays, bridge piers, flows around high-rise buildings, and a1111111111 a1111111111 a1111111111 a1111111111 a1111111111 quasi-periodic-I and quasi-periodic-II flow patterns. Moussaoui et al. [13] numerically examined the flow characteristics and heat transfer past three staggered square cylinders using the MRT-LBM. They studied the combined effect of Reynolds numbers (10 R e 100) and gap spacing (1 g 2) between the cylinders. They did not study the jet form between the cylinders and the presence of primary and secondary cylinder interaction frequencies for smaller and larger gap spacings. Alam and Zhou [14] experimentally examined the wake structure mechanism behind two identical side-by-side square cylinders at R e = 300 and 1 g 5. They discussed in detail the gap vortices, flow switch, merging of two streets into one and stability. Furthermore, they found single bluff body, narrow and wide streets, transition and the coupled-sheet flow patterns. Abbasi et al. [15] examined the effect of Reynolds number on flow past four square cylinders in an inline square configuration and discussed in detail the observed flow patterns. Chatterjee and Biswas [16] reported the chaotic and synchronous flow patterns at 1 g 5 for a fixed Reynolds number of 100. They argued that these flow patterns are the results of primary vortex shedding frequency and secondary cylinder interaction frequency for flow around two rows of staggered cylinders using Finite volume based computational fluid dynamics (CFD) solver using PISO algorithm. One can find the effect of Reynolds numbers and gap spacing for flow past three side-by-side square cylinders numerically [17][18][19]. These researchers discussed the effect of Reynolds numbers and gap spacing, and observed various flow patterns. Kumar et al. [20] numerically studied quasi-periodic, synchronous, and chaotic flow patterns at 0.3 g 12 for a fixed Reynolds number equals to 80.
The above literature shows that not much attention has been paid to flow past two rows of staggered square cylinders; in particular, the effect of Reynolds numbers has not been studied. This study aims to partially fill this gap in literature. Chatterjee and Biswas [16] have investigated the different flow patterns and there study is limited to a single Reynolds number and fixed the stream-wise gap spacing between the cylinders. The main goal of the present study is to systematically investigate the different flow patterns by changing the Reynolds number. We show the effect of Reynolds number from 20 to 140 with 1 g 6. The critical Reynolds number for the onset of vortex shedding frequency for the proposed problem is also discussed for the first time. In literature merging of jets [11,12,20] and variation in the wake size [13,14,17,18,19] has been considered separately whereas we discuss in detail the effect of jets between the cylinders on the wake interaction.
In Section 2 we provide problem description and numerical details. The grid independence, effect of the computational domain, and a validation are provided in Section 3 based on the available data of the single and two side-by-side square cylinders in literature. The dependence of flow patterns on R e and g is presented and discussed in Section 4 systematically. The main findings of this work are finally concluded in Section 5.

Problem description and numerical details
The purpose of this study is to simulate and analyze the flow characteristics of flow past two staggered rows of square cylinders Fig 1(A). The cylinders are fixed and identical in size (d) and shape. The spacing (s) between the cylinders is varied (g = s/d = 1, 2, 3, 4, 5 and 6) and the Reynolds numbers from 20 to 140. The upstream (C 11 , C 12 , C 13 , C 14 and C 15 ) and downstream (C 21 , C 22 , C 23 , C 24 , C 25 and C 26 ) rows have comprised five and six cylinders, respectively. The cylinders have been marked as C 21 through C 26 starting from the bottom of the computational domain. The stream-wise and the transverse directions with the origin of the coordinate system are x and y, respectively. The length of the computational domain L x = L u + d + s + d + L d (where L u = 8d is the upstream distance from the channel inlet and L d = 35d is the downstream distance from the downstream row of cylinders rear to exit of the domain) is taken as 45d. The height (L y ) of the computational domain is 11d. The velocity vectors are represented by u = (u, v). The simulations for the gap spacing 1, 2, 3, 4, 5 and 6 have been carried out using grids of sizes 441 × 921, 541 × 941, 641 × 961, 741 × 981, 861 × 1001, and 981 × 1021, respectively. An overall mesh at (R e , g) = (140, 4) is shown in Fig 1(B) as a representative case.
A uniform inflow velocity (u = U 1 ; v = 0) is prescribed at the inlet boundary of the channel. The convective boundary condition (@ t u + U 1 @ x u = 0) is applied at the exit of the computational domain [21]. The exit of the computational domain is far enough from the cylinders to ensure that there are no effects of the unconstrained movement of fluid. The periodic boundary conditions [22] have been applied on the bottom and top boundaries of the computational domain. No-slip boundary conditions (u = v = 0) are applied on the surfaces of the cylinders [21]. The forces on the surfaces of the cylinders are calculated using the momentum exchange method [23]. The simulation starts with the condition that the flow is initially at rest. The initial transient variations of physical parameters such as mean drag coefficient (C dmean ) and Strouhal number (S t = f s d/U 1 , where f s is the vortex shedding frequency) are not included in the obtained results. The S t is computed by performing Fast Fourier Transform (FFT) on the lift coefficient.
Two-dimensional lattice Boltzmann method (LBM) is used in this particular study. For an incompressible unsteady flow, the flow field can be described by the continuity equation and momentum equations, where u, ρ, p and t are the flow velocities, density, pressure and time, respectively. Applying the Chapman-Enskog expansion [24] with the restriction of low Mach number the above Eqs (1) and (2) can be obtained by the discrete lattice Boltzmann Eq (3) where f i (x, t) is the distribution function which indicates the position x of a particle at time t, Δt is the time step, f ðeqÞ i is the equilibrium distribution function and τ is the single relaxation time parameter which control the stability of the method. The equilibrium distribution function is defined by:  The kinematic viscosity ʋ can be obtained in the following way: The pressure, density and flow velocity can be obtained by The chosen range of Reynolds number fulfills the two-dimensionality requirements. In LBM, the movement of a large number of particles on a lattice is the basic proposition. In this study we have used D2Q9 (where D and Q denotes the dimensions and number of velocity particles, respectively) LBM model. Basically, LBM consist of two main steps: (i) streaming (left hand side of Eq (3)) and (ii) collision (right hand side of Eq (3)). The collisions take place between particles at each time step and their velocities change their directions but, the net mass and momentum are conserved. The boundary conditions are applied after the streaming step and the entire process is solved iteratively until the convergence is ensured by the follow- The LBM is found to be a second-order accurate method in space and conditionally stable depending on the single-relaxation-time parameter 'τ'. The most important advantage of this method is the ease of introducing the obstacle in the flow which makes this method fully suited for the present study. Here, the equilibrium distribution is found to be valid only for small Mach numbers. Also the density function is represented by particles that move one lattice length at every time step. Thus the speed of sound calculated from the diffusion velocity of particles therefore is always larger than the macroscopic velocity which results in small Mach number restriction. The lattice Boltzmann method actually belongs to a class of the pseudocompressible solvers of the Navier-Stokes equations for incompressible fluid flow. In order to correctly simulate the incompressible fluid flow, it must be ensure that the Mach number, M a = U 1 /C s << 1. Furthermore, low Mach number requires that the uniform inflow velocity U 1 be sufficiently small.
The present computational code is originally written and edited using the Intel Fortran platform with 64-bits. All the computations were performed on an Intel i7-2600, four core 3.2GHz, 4GB DDR3 memory, 320GB hard disk computer, with a Window 7 (64-bits) system. The present simulation adopts the uniform and square grid. More details can be seen in review article of Chen and Doolan [25]; books of Sukop and Thorne [22]; Succi [26] and Kruger et al. [27], and research articles [28,29].

Grid independence, effect of computational domain and code validation study
The basic validation of this present computational code, basics of lattice Boltzmann method, grid independence study and effect of computational domain has already been carried out for flow past a rectangular cylinder [30], and flow past more than two side-by-side cylinders [18,19,28,29]. In Fig 2(a) and 2(b) and Tables 1 and 2 the C dmean1 and C dmean2 are the mean drag coefficients of the lower and upper cylinders, respectively, in two side-by-side square cylinders arranged from bottom to top. Similarly, S t1 represents the Strouhal number of lower, and S t2 represents the Strouhal number of upper cylinders. The drag and lift coefficients, C d and C l , are obtained by normalizing the forces by 2F d /ρU 2 1 and 2F l /ρU 2 1 (ρ is the fluid density, F d is the force along stream-wise direction and F l is the force along cross-stream-wise direction).
A grid independence study was performed by comparing the calculated flow characteristics using 10, 20 and 40 points to discretize the cylinder side for flow past two side-by-side square      cylinders. A maximum difference of 1.3%, 1.3%, 0.6% and 0.6% on C dmean1 , C dmean2 , S t1 and S t2 , respectively, was found between the last two cases. In this study, we adopt a resolution of 20 points for square cylinder. In addition, we choose 20 points, keeping in mind the constraints of long running simulations and computational resources. Kumar et al. [20] found that 16-20 points cylinder resolution is enough to obtain grid independent results. We also found that L u = 8d, L d = 35d and L y = 11d is the best choice in terms of computational length see Table 1.
In addition, for the sake of the present numerical study, we have carried out a separate validation study for flow past a single square cylinder see Fig 3(a), 3(b), 3(c) and 3(d) and Table 2 and two side-by-side square cylinders (see Table 3) at (R e , g) = (73, 4). It is to be noted that in case of single square cylinder in a channel we adopted the no-slip [21] boundary conditions on the top and bottom walls of the channel. The same boundary conditions and channel length adopted as discussed already for the proposed problem in Section 2. The vortex pattern in the wake behind the cylinder shown in Fig 3(a) composes of negative and positive alternate shed vortices, i.e. the well-known Karman vortex street at R e = 140. The velocity profile is also presented in Fig 3(b). It is observed that the drag and lift coefficients are settles to a periodic behavior Fig 3(c). Furthermore, no other peaks observed in the power spectra analysis of lift coefficient for R e = 140 in Fig 3(d). In this paper in spectra plots 'PSE' stands for power spectrum energy. These findings are in good agreement with those observed by Davis and Moore [31]. Similar characteristics observed for R e = 75, 100 and 120 (not shown). The force statistics data of different Reynolds numbers are given in Table 2 with available data in the literature [32,33]. These findings further gives us the confidence that the present code recovered correctly all the silent features of an isolated cylinder in a cross-flow. The comparison of present results and those of Agrawal et al. [11] are given in Table 3. The minor discrepancy in S t found could be due to the difference in the cylinder size and computational domain, etc. used.

Results and discussions
To the authors' knowledge, there are no experimental measurements and numerical studies of the flow patterns for the problem under consideration in the open literature to compare with the numerical work of this study. We will mostly compare our results with flow past a single, two, three, and row of square cylinders available in the literature. Results are presented for two-dimensional computations for 75 R e 140 and 1 g 6 using the single-relaxationtime lattice Boltzmann method. In this study C 12 , C 13 , C 23 , and C 24 of two consecutive cylinders (arbitrarily chosen) for time history analysis of drag and lift coefficients, and power spectra of lift coefficients to further verified the shedding pattern of observed flow patterns. In vorticity contours visualization plots the solid and dashed lines represent the positive and negative shed vortices generated from the lower and upper corners of the cylinders. We used the solid and dotted lines in drag and lift coefficients plots. The solid line used in power spectra graphs for chosen cylinders in this study.
One typical example of onset of vortex shedding of flow from steady to unsteady is shown in Fig 4(a) and 4(b) at g = 3 for two different Reynolds numbers. Our study found that the maximum value of critical Reynolds number (R ecr ) is 41.5 at g = 6, and reduces with a reduction in gap spacing see Fig 5. It is to be noted that to find the critical Reynolds number for each gap spacing we start from R e = 20 and increase it with the increment of 0.5 until we not find the critical value. This is actually due to increase in the vorticity production from the boundary layers with decreasing the separation between the cylinders as a result of gap induced accelerated effects. With a reduction in gap spacing the vorticity generation increased is the main cause of the reduction of critical Reynolds numbers with increasing the gap spacing [12]. Therefore, the lower limit of Reynolds number for further analysis is set at 75. Saha et al. [34]  Flow patterns and aerodynamic characteristics around staggered rows of cylinders observed that the flow loses its two-dimensionality between 150 R e 175, therefore the upper limit of R e in this study has been decided as 140.
The instantaneous vorticity contours visualization plot at (R e , g) = (75, 6), (140, 6), (100, 5) and (120, 4) are presented in Fig 6. For relatively large gap spacing, g = 6, 5 and 4 the shed vortices are clearly seen from both rows of the cylinders. For (R e , g) = (100, 5) and (120, 4), the merging of vortices observed at the downstream region; however, for (R e , g) = (75, 6) and (140, 6) the shed vortices throughout the computational domain almost remain distinct. Fig 6  (a)-6(d) represents the change from the anti-phase to in-phase mode. At the beginning, the streets behind the cylinders are anti-phase, characterized by the same lateral spectrum between the shed vortices and same width. As a result, the vortex shedding frequencies from the cylinder differ. It is observed that anti-phase vortex shedding is predominant for larger gap spacing between the cylinders as compared to in-phase vortex shedding behind two consecutive cylinders. Williamson [35] found experimentally the anti-phase vortex shedding flow pattern for Flow patterns and aerodynamic characteristics around staggered rows of cylinders two circular cylinders at (R e , g) = (100, 3). Chatterjee and Biswas [16] also had similar numerical observations. The time history analysis of drag and lift coefficients are given in Fig 7(a)-7 (p). By changing the Reynolds number one can clearly observe the significant variation in the signal. The lift coefficient signals in these figures clearly indicate that the secondary cylinder interaction frequency almost vanishes. The shedding pattern behind the cylinders shows some in-phase and some anti-phase characteristics. Qualitatively, the flow behavior is in good agreement with the results obtained numerically by Chatterjee and Biswas [16] for flow past two staggered rows of square cylinders. This can be further revealed from the power spectra analysis of lift coefficients presented in Fig 8(a)-8(p), where S ts = 0.3047 and 0.3006 for (R e , g) = (75, 6) of C 23 and C 24 , respectively (where subscript 's' denotes secondary). Similarly, one can see S ts for C 23 and C 24 , respectively; in these figures. The power spectra analysis of lift coefficients clearly indicates the dominancy of primary vortex shedding frequency (S tp = 0.1518 and 0.1518 for C 12 and C 13 , respectively; for (R e , g) = (75, 6); where subscript 'p' denotes primary). This is interesting that the secondary cylinder interaction frequency disappears at the larger Reynolds number for g = 6, 5 and 4. The concept of secondary cylinder interaction frequency was proposed by Kumar et al. [20]. They proposed that the existence of Flow patterns and aerodynamic characteristics around staggered rows of cylinders secondary cylinder interaction frequency occurs due to narrowing and widening of the wake behind one cylinder. Chatterjee and Biswas [16] in there simulation for (R e , g) = (100, 5) obtained S tp = 0.187, which fairly agrees with our present computed results although, the Reynolds numbers are somewhat different. The downstream row of cylinders is shedding vortices as an isolated cylinder and no sign of the frequency of oscillation of the upstream row of cylinders are noticeable in the spectra analysis of lift coefficients. Our results are in good agreement with those of the Kumar et al. [20] for flow past a row of square cylinders. . On the basis of these characteristics, we called it quasi-periodic-I flow pattern. Sewatkar et al. [12] had similar observations at (R e , g) = (100, 3) for flow past a single row of square cylinders.
The instantaneous vorticity contours visualization plot at g = 2 for R e = 75 and 100 are shown in Fig 12(a) and 12(b). For R e = 100, the shed vortices are seen approximately up to x/d % 10 and then coalesces with each other shed vortices further downstream of the computational domain. One can clearly see the phase difference of shedding from an adjoining cylinder of second row. It is observed that the shedding phase difference of certain adjoining cylinders is maintained for larger time. In case of two side-by-side circular cylinders for R e = 100 and 1 g 5, Williamson [35] observed that the changeover between in-phase and anti-phase vortex shedding mode lasts for a large number of vortex shedding cycles. Furthermore, at g = 2 the time history analysis of drag and lift coefficients for R e = 75 in Fig 13(a)-13 (d) reveals that the primary and secondary cycles are not showing constant behavior for different time periods. Similar drag and lift coefficients characteristics are observed in Fig 13(e)-13 (h) for R e = 100, respectively. The power spectra analysis of lift coefficients at g = 2 for different Reynolds number indicates S tp = 0.0.0149, 0.0089, 0.2160 and 0.2113, and S ts = 0.0208, 0.0327, 0.2101 and 0.1994 for C 12 , C 13 , C 23 and C 24 , respectively; for R e = 100. On the basis of above characteristic, the flow is so-called quasi-periodic-II flow pattern in Fig 14(a)-14(h). In quasiperiodic-II flow pattern both the primary and secondary cycles have variable time period.
At g = 1, the instantaneous vorticity contour plots for different Reynolds numbers are shown in Fig 15(a) and 15(b). The shed vortices coalesce together just near downstream of the second row of cylinders for R e = 100. This shedding coalesce occurs at small downstream of the second row of cylinders and some of the wakes appear wider and narrower. As a result, the strong merging and distortion of shed vortices clearly can be seen as they move downstream throughout the entire computational domain and the flow completely behaves chaotic. The similar characteristics observed for R e = 140 in Fig 15(b). The time history analysis of drag and lift coefficients for different Reynolds numbers confirms that the flow is chaotic in nature which is presented in Fig 16(a)-16(h). The lift forces on downstream row of cylinders in these figures are very irregular and anti-phase with each other (only occasionally). It is observed for Fig 13. Time-history analysis of the drag and lift coefficients (a-d) (R e , g) = (75, 2) and (e-h) (R e , g) = small gap spacing that the flow interference between successive cylinders is too strong. Furthermore, the jets are merging immediately behind the cylinders and also interact with the shed vortices. These observations are also verified from the time history analysis of the drag and lift coefficients, and the power spectra analysis of lift coefficients (Fig 17(a)-17(h)).The power spectra analysis of lift coefficients further reveals the behavior of chaotic flows as a result of broad and continuous spectrum. It is concluded that at lower gap spacing the shed vortices interact so strongly that the Reynolds number does not affect the basic characteristics of the flow from R e = 75 to 140. We further observed for chaotic flow pattern that there is no relation between shed vortices from different cylinders and showing haphazard motion of shed vortices from the cylinders. The secondary cylinder interaction frequency completely dominates the primary vortex shedding frequency. The chaotic flow pattern was further confirmed from the phase portrait plots shown in Fig 18(a)-18(h). From these plots one can clearly see that the dynamics of the force coefficients of two adjoining cylinders progressively decouples, as suggested by Ravoux et al. [36]. The chaotic nature of flow is shown in the portrait diagram in these figures of drag versus lift coefficients for (R e , g) = (140, 1). The shed vortices behind the Flow patterns and aerodynamic characteristics around staggered rows of cylinders cylinders are strongly decoupled. However, the strong jet flow between the cylinders causes a chaotic state.
In synchronous flow pattern, the lift forces are almost symmetric, resulting in zero mean lift force as observed in portrait plots. The periodic and repulsive forces between the downstream row of cylinders have higher amplitude than those between the upstream row of cylinders. The irregularity of the lift coefficient on downstream row of cylinders is due to the jet flow between the cylinders. The trajectories in portrait plots are remains in figure-eight shape for last four figures.
The relationship between the statistical characteristics and vortex patterns is another interesting issue and will be discussed in detail in this section. In coming figures the aerodynamic characteristics are presented in mean drag coefficient, Strouhal number, root-mean-square values of drag and lift forces. For comparison, the results of single isolated cylinder are also given in these figures. It is found that the C dmean depends both on the R e and g. The C dmean of upstream and downstream row of cylinders are quite different. The C dmean for the very small gap spacing (g = 1 and 2) is larger compared to the other gap spacing between the cylinders. C dmean for isolated cylinder is smaller than the flow past two staggered rows of cylinders configuration, except the case of (R e , g) = (75-100, 1), where C dmean of the cylinders (C 25 , C 26 ) is lower than the value of single square cylinder. As observed in Fig 19(a)-19(h), the wake response transition behind the downstream row of cylinders leads to significant change in C dmean . For (R e , g) = (75-100, 2), a sudden jump in C dmean of C 25 and C 26 is due to changeover of the chaotic flow pattern to quasi-periodic-II flow pattern. As the gap spacing between the cylinders is increased to g = 6 and 5, the synchronized vortex shedding behind the cylinders dominates the wakes, then the C dmean of downstream row of cylinders is almost similar to that for flow past a single square cylinder. However, it is not the case for upstream row of cylinders.  Flow patterns and aerodynamic characteristics around staggered rows of cylinders and very close to isolated cylinder value, for example, S t = 0.1392, 0.1498, 0.1542, and 0.1565 at R e = 75, 100, 120, and140, respectively. Sudden jumps in vortex shedding frequency are observed for quasi-periodic-II and chaotic flow patterns. We have not observed Strouhal value for those cylinders which shows steady behavior. The Strouhal number values for individual cylinder are smaller than that of a single cylinder for downstream row of cylinders at g = 1 for all chosen Reynolds numbers. On the other hand, in the case of upstream row of cylinders it is smaller at g = 2 for R e = 75; g = 1 and 2 for R e = 100 and 120; and g = 2 and 3 for R e = 140. The discontinuities in Strouhal number occur at the Reynolds number between any two flow patterns. When the wake of the two rows of staggered square cylinders transits from chaotic flow pattern to quasi-periodic-I flow pattern at (R e , g) = (75, 1-2) and (100, 1-2); and from quasi-periodic-I flow pattern to quasi-periodic-II flow pattern at (R e , g) = (75, 2-3) and phase mode. Dramatic drop happens similar to that of quasi-periodic-II flow pattern at (R e , g) = (100, 3) and quasi-periodic-II flow pattern at (R e , g) = (100, 2).
The root-mean-square values of lift coefficients are given in Fig 22(a)-22(h). The single square cylinder values are also given in these figures for comparison. The very small C lrms for g = 6, 5 and 4 indicates that the vortex shedding is very weak behind the cylinders for synchronous flow pattern. Furthermore, the C lrms of downstream rows of cylinders are close to isolated cylinder with the increase of gap spacing. At all these relatively large gap spacing (g = 6, 5 and 4), we found that as the Reynolds number changes, transition in flow patterns from quasi-periodic-I to synchronous occurs (Fig 23).
The flows patterns map have been demarcated based on the characteristics observed in terms of vorticity contours visualization; time history analysis of the drag and lift coefficients; and power spectra analysis of the lift coefficients. The flow patterns map as a function of R e and g is given in Fig 23. As we observed synchronous, quasiperiodic-I, quasiperiodic-II and chaotic flow patterns for 75 R e 140 and 1 g 6. These flow patterns were also observed by Sewatkar et al. [12] and Kumar et al. [20] for flow past a single row of square cylinders. The synchronous flow exists at higher Reynolds number at g = 4 and larger gap spacing (g = 5, 6) for all chosen Reynolds numbers. The quasi-periodic-I flow pattern was observed at g = 4 and Flow patterns and aerodynamic characteristics around staggered rows of cylinders 3 with R e 100 and R e 140, respectively. Quasi-periodic-II flow pattern was found at g = 2 for Reynolds numbers 75 and 100. The figure also depicts the chaotic flow pattern for all chosen Reynolds numbers at smaller gap spacing (g = 1) and at (R e , g) = (120, 2) and (140, 2). The same signed vortices and its strength were found behind the second row of cylinders and approximately equal to that observed for flow past a single row of square cylinders in the synchronous and quasi-periodic-II flow patterns by Sewatkar et al. [12]. As a result, the shed vortices move independently behind the second row of cylinders in stream-wise direction. In Quasi-periodic-I and chaotic flow patterns significant merging of shed vortices and transverse movement were observed behind the second row of cylinders throughout the computational domain. We found that both the Reynolds number and gap spacing significantly affects the flow characteristics, with the effect of gap spacing greater than that of the Reynolds number.

Conclusions
The numerical results of two-dimensional flow across two rows of identical square cylinders in staggered arrangement are reported in this study using a lattice Boltzmann method. The main goal of the study is to systematically understand the Reynolds number effect on steady to unsteady transition and to further investigate the different flow patterns, at certain spacing. Furthermore, we study the existence of secondary cylinder interaction frequency in the time Flow patterns and aerodynamic characteristics around staggered rows of cylinders series of lifts coefficients. It is found that the critical Reynolds number, at which onset of vortex shedding occurs, increases with an increase in gap spacing between the cylinders. It is also observed that the Reynolds number have substantial effect on the flow characteristics especially at g = 4 and 2. For g = 6 and 5, the secondary cylinder interaction frequency disappears at large Reynolds numbers and the flow completely dominates by the primary vortex shedding frequency. This confirms that at such large gap spacing the wakes interaction behind the cylinders is weak and it reduces further with an increase in Reynolds number. For g = 1 and 2 the secondary cylinder interaction frequency predominates the flow and as a result strong wakes interaction are observed behind the cylinders at all Reynolds numbers.
The wake interaction mechanism is also discussed in detail. At larger gap spacing (g = 6 and 5), we found weak wakes interaction and throughout the entire computational domain between the cylinders a continuous jet forms. At relatively small gap spacing (g = 4 and 3), the wakes interaction behind the cylinders strongly depend on the Reynolds number. The wakes interaction occurs due to jets coming out of the spacing between the cylinders. In quasi-periodic-I flow pattern, the secondary cylinder interaction frequency depends on the spacing between the cylinders. In quasi-periodic-II (g = 2) flow pattern, the jets effects is stronger. In chaotic flow pattern (g = 1) the jets effect is too strong and the wakes are immediately broken behind the cylinders. The primary vortex shedding frequency is strongly affected by the secondary cylinder interaction frequency at smaller gap spacing.
Finally, in this work we clearly brought out the wakes interaction behind the cylinders and the importance of jets forms between the cylinders which were observed experimentally for higher Reynolds numbers. We also discussed in detail that the physical parameters also affected due to changeover of flow patterns. Flow around a two rows of staggered square cylinders at higher Reynolds numbers will be studied in near future using the three-dimensional lattice Boltzmann simulations.

Author Contributions
Formal analysis: Shamsul Islam.