Computational analysis of naturally oscillating tandem square and circular bluff bodies: a GPU based immersed boundary – lattice Boltzmann approach

ABSTRACT In this study, a fluid-structure interaction solver is employed to explore the flow physics of tandem oscillating square bluff bodies by varying the corner radius. Based on previous research, dynamic behavior, and fluid forces only square and circular bluff bodies are chosen for a supplementary investigation. This study investigates the effects of the streamwise displacement, , of two elastically mounted square and circular bluff bodies in a tandem arrangement. Fluid flows at a Reynolds number of 100 and a fixed mass ratio are used for the upstream and downstream bluff bodies. As the most complex oscillations and wake flows occur adjacent to the approximate natural frequency ratio, thus a constant natural frequency of is defined. To achieve excellent parallel performance, an immersed boundary lattice Boltzmann method code coupled with a structural equation is developed and accelerated using a graphical processing unit. The results for flow-induced vibrations suggest that the flow response characteristics of the tandem oscillating bluff bodies are strong functions of the streamwise displacement. Four different flow regimes of the bluff body arrangements are identified on the basis of response amplitude. Both the square and circular downstream bluff bodies exhibit two peak transverse amplitudes in the existing scrutinized spacing ratios. In addition, the critical spacing ratio values are dissimilar for the two bluff bodies. Moreover, at the end of the spacing ratio ( ), the wake effects for the downstream square body are amplified, and it oscillates near a single oscillating circular bluff body. Finally, the flow physics, as well as the hydrodynamic force coefficients of the tandem oscillating bluff bodies are also reviewed.


Introduction
The interactive fluid-structure phenomena of bluff body wake and flow separations have practical importance and vital significance in considerable scientific disciplines, such as hydrodynamics, aerodynamics, thermodynamics, biology, civil engineering, and so on. The wake produced by a nonlinear interaction between a fluid and a flexible mounting structure moves the structure under the fluid force as vortex shedding frequency matches the structure frequency. This phenomenon is called flowinduced vibration (FIV).
The field of FIVs significantly widened after the disaster of the Tacoma Narrows Bridge in 1940. Practically, vibration-related problems cause wide-ranging and expensive damage to structures, reduce their lifetimes by fatigue, and pose risk to structures. Païdoussis (2006) summarized some of the vibration-related problems during nuclear power plant operation. Currently, there are various technological concepts, and one of them is harvesting energy from FIVs, which is a promising approach. With the practical approach of obtaining green CONTACT Hojin Ha hojinha@kangwon.ac.kr; Ehsan Adeeb ehsan_adeeb@hotmail.com energy, Bernitsas et al. (2008) proposed a device design (VIVACE) that converts spring-mounted circular bluff body oscillations into electrical energy under the action of a water current.
It is well known that as fluid flows around a bluff body, tremendous changes occur in the flow and wake regions with the variation in the Reynolds number. Moreover, flow separation plays a significant role in bluff bodies, and the point of separation for a square body is fixed by its leading and trailing edges, whereas it is not fixed for a circular body (Sahu & Patnaik, 2011;Sharma & Eswaran, 2004). Concerning the streamwise displacement and the Reynolds number, numerous studies have found three regimes for tandem stationary square and circular bluff bodies: single bluff body regime, reattachment regime, and co-shedding regime. The regime streamwise displacements are different for tandem stationary square and circular bluff bodies. Moreover, based on the recent studies by Sumer and Fredsøe (2006) and Adeeb et al. (2018b), drag inversion/critical spacing phenomena are specifically sensitive to the Reynolds number, bluff body shape, and streamwise displacements (3.0D-5.0D) of bluff bodies.
If one or more bluff bodies are attached to a spring and a damper and allowed to move in a longitudinal direction, they undergo FIVs, wake-induced vibrations (WIVs), and wake-induced galloping (which is not covered in this study). Lock-in, a well-known phenomenon, yields large amplitudes for bluff bodies when the bluff body oscillation frequency matches the vortex shedding frequency (Williamson & Roshko, 1988). Both numerical and experimental techniques have been employed to study FIVs of single and tandem bluff bodies. Miran and Sohn (2017) employed a finite-volume code to study the effect of the corner radius of a square bluff body for a specific Reynolds number of 100. Several computations were conducted by varying the frequency ratio. In their study, the minimum response was detected for a single square bluff body, and its counterpart single circular bluff body exhibited the maximum response amplitude. Adeeb and Sohn (2021) also studied FIV phenomena for a square bluff body by varying the corner radius. The obtained results presented trends similar to the above. Specifically, the lock-in region for the square bluff body was very narrow; however, it expanded as the corner radius changed. The widest region is observed for a circular bluff body, and all other transverse amplitude regions lie between the two bluff body shapes. Hong et al. (2021) used the ghost-cell immersed boundary method to simulate moving structures with finite or zero thickness. To obtain the stable solution for the large time step size of the dynamic thin elastic structure, they employed the composite implicit time integration method. The present developed method is capable of capturing the flow physics of several canonical examples with mass conservation. Furthermore, various other numerical and experimental studies have been conducted to understand FIV phenomena. Meneghini and Bearman (1999), Bokaian and Geoola (1984), Akbari and Price (1997), Sumer and Fredsøe (2006), Triantafyllou et al. (2016), and Anagnostopoulos (1997) proposed that the response amplitudes and forces on bluff bodies increase for a range of frequencies and the maximum values are shown in the lock-in region. Zhou et al. (2013) performed numerical simulations to determine the response amplitude of a square bluff body and presented the effect of the angle of incidence. The minimum response amplitude and its corresponding range of lock-in were observed at zero angles of incidence. Hines et al. (2009) numerically investigate the vortexinduced vibration of a square cylinder by introducing the localized mesh deformation method to the STREAM code. They witnessed Hysteresis phenomena within the lock-in range related to phase difference. In addition to that, they simulate an autonomous vibration motion model, whose results are evidence that a coupled response exit between the square cylinder and flow field. Huera-Huarte and Beeraman (2011) investigated the dynamic responses of tandem circular bluff bodies in a water channel with low mass ratios (1.45 and 1.80) and in close proximity (2.0 < L x /D < 4.0). They explained that the dynamic response amplitude of the downstream bluff body was composed of either FIVs, WIVs, or a combination of FIVs and WIVs. Mahir and Rockwell (1996) conducted experiments on a small-Reynolds number flow (Re = 160) with tandem oscillating circular bluff bodies. The focus was on the effects of the bluff body response amplitudes and the phase angles for a wide range of frequencies. It was concluded that the lock-in regime for the tandem bluff bodies was wider than that for a single bluff body. Moreover, as the separation distance reduced, the range of frequencies at which lock-in occurs also increased.
In addition to experimental research, numerous numerical approaches have been applied to study FIV phenomena for tandem square and circular bluff bodies. Mital et al. (1997) conducted pioneering research in utilizing computational resources to examine FIV phenomena for tandem circular bluff bodies. Two specific streamwise distances (2.5D & 5.5D) were chosen for their study, and the Reynolds numbers were set as 100 and 1000. They reported that the results depended significantly on the streamwise position and the Reynolds number. Furthermore, as the wake of the upstream bluff body approached the downstream bluff body, it started to exhibit a significantly large-amplitude owing to a wakeinduced flutter. Papaioannou et al. (2008) employed the spectral element method and set the Reynolds number for the fluid as 160. Three different streamwise positions were selected in the range of 2.5 < L x /D < 5.0. Based on the spacing, a shift in the reduced velocity was observed for the upstream bluff body. Furthermore, an upper branch and hysteresis phenomena were only observed for the single oscillating bluff body; as the second bluff body was introduced downstream, these effects did not occur for it. Recently, E. Wang et al. (2017) numerically investigated the FIVs of two identical oscillating circular bluff bodies using a two-way fluid-structure interaction (FSI) method at a Reynolds number of 500. For the downstream bluff body, two response regimes as classified by other studies were identified, whereas the third regime was missing. One possible reason for this was that they fixed both ends of the downstream bluff body.
To avoid geometrical complexities, researchers and designers prefer streamlined shapes in their engineering practices, such as square and circular bluff bodies.
In contrast to dual circular bluff bodies, two square bluff bodies have also generated research interest, the literature on this topic is limited and lesser than that on circular bluff bodies. Kumar and Sen (2021) numerically studied the wake interactions of tandem square bluff bodies for a fixed spacing ratio of 5.0D by varying the reduced velocity and employing a stabilized space-time finite-element solver. Their simulations showed that the frequency characteristics of the two bluff bodies were identical; however, the upstream bluff body separately considerably compared with its downstream counterpart. Han et al. (2018) numerically examined tandem oscillating square bluff bodies (upstream static) in both transverse and streamwise directions for a wide range of reduced velocities. They identified six vortex modes 2S, 2S * , 2P, P + S, 2T, and steady wake modes for all ranges of the Reynolds number and the reduced velocities.
In recent years, flow control has been a major focus in research, leading to substantial benefits in industrial applications. In this research, a passive method is applied, where the corner radius of a square bluff body is altered to eventually achieve a circular body (Adeeb et al., 2018a). Furthermore, the current study is restricted to square and circular bluff bodies.

Numerical methodology
Various numerical techniques are available for solving computational fluid dynamics problems. In this research mesoscopic approach is applied: where the behavior of a collection of particles is represented by a distribution function. The kinetic theory is the main support of the mesoscopic method, on which the lattice Boltzmann method (LBM) is based.

Lattice Boltzmann method (LBM)
In perspective of the current study, the flow field is governed by the discrete lattice Boltzmann equation (LBE) with a single relaxation collision operator (Bhatnagar et al., 1954). The discretized lattice Boltzmann equation is solved at each time step, which is written as In Equation (1), t is the time, c is the particle velocity, and f is the particle distribution function, f a (x, t) denotes the density distribution function at time step t, and at the next time step, t + t, the distribution function moves to the adjacent neighboring node, x + c a , with discrete velocity c a . For 2D problems, most commonly used model is D2Q9, in which the lattice is spread in 2D with nine velocity vectors.
In the lattice Boltzmann framework, macroscopic variables, such as density and velocity, at each lattice site are defined using the local density distribution function, as expressed in Equation (2). In addition, the fluid pressure is determined using the equation of the state. The implementation of a no-slip boundary condition in the LBM is complex, and its significance should be considered carefully. Boundary conditions influence the entire computational domain, although they are applied to a very small portion of it. In the LBM, no-slip boundary conditions are applied indirectly. This requires an evolution of distribution functions that satisfy the boundary conditions. The simplest approach for representing any obstacle in the LBM is the staircase approximation (bounce-back boundary condition), as shown in Figure 1. Moreover, accuracy can be achieved by increasing the mesh size; however, this is at the cost of large computational resources and a long computational time.

Immersed boundary method (IBM)
The IBM (immersed boundary method) method is used in the LBM to overcome the problems of a high skew mesh, a negative cell volume, and non-orthogonality, specifically in body-fitted grids. The main concept of the IBM is that body markers can be set off-lattice, which affects the fluid by a force field, which is computed using an interpolation stencil. Therefore, the fluid mesh does not need to confine the solid obstacle, and the two domains can operate independently. Peskin (1972) was the first to propose the IBM model to simulate blood flow patterns around heart valves. As shown in Figure 2(a), in the IBM, two sets of points store information: Lagrangian points, which represent complex geometric shapes, and an Eulerian grid, on which the governing equation of the flow field is solved. The interaction between these two sets of points is computed using various interpolation schemes. Because the Lagrangian points are in motion, the IBM only tracks the moving body data points, with the Eulerian grid remaining stationary, which saves significant computational time.
In this study, for IBM, an approach similar to that used by Adeeb et al. (2018b) is used. To mimic a no-slip boundary condition, a velocity-correction-based IBM (X. Wang et al., 2014) is coupled with the LBM. In the formulation, the velocity is modified at the neighboring points of the immersed boundary. The corrected velocity  field is obtained by Lagrangian interpolation to satisfy a no-slip boundary condition at the boundary intersection points (px, py) ( Figure 2(b)), as expressed in Equation (3) where k denotes the velocity corrected points. a and b at the horizontal mesh line, c and d at the vertical mesh line, i denotes the points in the x-direction (p x , a, a 1 , b, and b 1 ), and j denotes the points along the y-direction (p y , c, c 1 , d, and d 1 ). The mean value is used if the velocity is modified more than once at a mesh point. Subsequently, the corrected velocity is used to determine the body force as follows (Equation (4))

Structural model
The structural dynamic system of a rigid body subjected to aerodynamic or hydrodynamic forces is represented by a mass-spring-damper system. A rigid body immersed in water/air functions as a mass-spring-damper system and mathematically can be represented in the transverse direction as in Equation (5) m ÿ + cẏ + ky = F where m is the structure mass, F is the force on the bluff body, k is the stiffness, and c is the damping coefficient. Here, [ÿ,ẏ, y] are the acceleration, velocity, and position of the bluff body, respectively. In the nondimensional form, the motion equations can be written as for the cross-flow direction (transverse) in Equation (6).
Above, F L is the force applied to the bluff body by the fluid in the transverse direction. The nondimensional parameters are defined in Equation (7) as where ζ is the damping ratio, ρ is the density of the fluid, m * is the mass ratio, U * is the reduced velocity, U ∞ is the free-stream velocity, f n is the natural frequency of the bluff body, f s is the vortex shedding frequency of the stationary bluff body, D is the diameter of the bluff body, and F r is the frequency ratio. The integration of the forces in Equation (8) is the resultant aerodynamic or hydrodynamic force exerted on the bluff body. The forces in the x-and y-directions can be obtained as where drag coefficient C D and lift coefficient C L are expressed as in Equation (9) After computing the forces at each time step, Equation (6) is solved using the fourth-order Runge-Kutta (RK-4) method. The time-averaged drag force and the rootmean-square lift coefficient values are calculated using Equation (10) The response of an oscillating bluff body to the timeaveraged transverse amplitude is defined as In Equation (11) y max and y min are the transverse maximum and minimum displacements, respectively.

Parallel computing
To further enhance this study, the code is extended for GPU (Graphics Processing Unit) parallel computing. One of the most significant aspects of compute unified device architecture (CUDA) is its execution model. Overall GPU execution is summarized as follows. First, the host launches a kernel in a device with the information of the grid and block size. Subsequently, block threads are assigned to the device multiprocessor, which executes all threads into wraps.

GPU algorithm
Hence, in the developed source code, threads are assigned such that each thread updates one lattice node and executes the same kernel in the computational grid.
Step 1 Space allocation: The first step is to allocate the memories on the device and the host. In the device memory, the thread ID is allocated as follows: where NX, NY are a total number of lattice sites in x and y directions, respectively, k is the number of linkages in the lattice model, and i and j are the thread unique IDs in the x and y directions, respectively.
Step 2 Initialization: All physical parameters are converted into lattice units, subsequently, the initialization kernel is called, which specifies the value according to the equilibrium distribution function.
Step 3 Macroscopic properties computation: This kernel calculates all macroscopic properties on each lattice node.
Step 4 Immersed boundary -The immersed boundary kernel computes the value of the force on each boundary point using the above-mentioned method. In addition, Equation (6) is used to update the rigid bluff body motion using the RK-4 model.
Step 5 Boundary function: This kernel consists of different boundary conditions (such as inlet, outlet, and periodic). In addition, each boundary condition is executed individually.
Step 6 Equilibrium: This kernel calculates the equilibrium value for each lattice node.
Step 7 Collide: This kernel executes the collision model (Single Relaxation Time, SRT) between molecules.
Step 8 Streaming: This kernel computes the standard advection to all lattice nodes depending on the lattice model. In streaming (the process of advection), the density distributions are moved to the adjacent neighboring lattice site depending on the lattice model. The streaming process is as follows (Equation (12)) Step 9 Step 3 is adopted until convergence. In the developed source code, after 1000 iterations, all data are transferred back to the host from the device.  In this research, computations are performed on an Intel core (i7-6700K, 4.00 GHz) CPU with 16 GB RAM. Three different types of GPU cards are used: GTX 970, GTX 1070 Ti, TITAN V and shown in Figure 3. The operating system is Windows 7, and the CUDA 8.0 toolkit is used.

Problem description
A schematic of the computational domain for two tandem square and circular bluff bodies of equal diameter D is shown in Figure 4. Air is selected as the working fluid, and all computations are performed for a fixed Reynolds number of 100. The inlet (boundary condition u = U ∞ , v = 0) is located at a distance of 10D from the upstream bluff body, and the pressure outlet (p = 0 (boundary condition)) is located at 30D from the downstream bluff body. Both lateral boundaries are at 12.5D from the center of the bluff bodies where the slip boundary condition is applied. Using the IBM, a no-slip boundary condition is applied on the surfaces of the bluff bodies. Numerous tandem configurations are investigated for both tandem bluff bodies (square and circular) by varying the centerto-center distance between them as 1.05 ≤ L x /D ≤ 12.0.
The structural parameters are chosen, such as a fixed mass ratio of m * = 10, and to promote high oscillations, a zero damping parameter (ζ ) is selected. Moreover, all simulations are conducted for a fixed frequency ratio F r = 1.0.

Parallel performance
An hybrid CPU/GPU approach is used in which the computationally extensive parts of the simulation are implemented on the GPU, whereas the remaining part is performed on the CPU. The performance directly depends on the capability of the GPU card, complexity in the LBM model, LBM velocity model, and some other parameters that are harder to control. The GPU performance is evaluated in terms of, the million-lattice site update per second (MLUPS). In a 2D simulation, the MLUPS is simply calculated from Equation (13) where NX * NY represents the size of the lattice site update and T is the time (in seconds). To assess the 2D LBM code performance, the GPU code is compared with a single-thread CPU code written in CUDA and C++ respectively. For the first case, the performance of a liddriven cavity flow problem for various grid sizes is analyzed. The GPU performance increases with the increase in the mesh size and ranges between 17 and 609 times faster than the CPU (improve with the increase of fluid data point density). The second test case corresponds to a flow over a circular bluff body at the Reynolds number of 100using the LBM. However, for the flow over a circular bluff body, additional boundary conditions need to be solved while executing the LBM algorithm. With the increase in complexity, the GPU performance decreases as the reaming threads become idle. The maximum performance is achieved from the TITAN V card, which is 145 times faster than a CPU code. For the IB-LBM case, the GPU performance is evaluated again on a flow over a circular bluff body at a Reynolds number of 100. Because of the hybrid approach, the IBM part is solved on the CPU. Table 1 summarizes the results of the hybrid approach. The increasing complexity of the LBM causes a greater delay in the execution time (averaged over 500 iterations). Furthermore, TITAN V outperforms all other cards easily, whereas it loses approximately 60% performance compared to the GPU-LBM code owing to branch divergence; however, it is still 90times faster than the CPU code.

Code validation
To validate the proposed FSI scheme, the simulated test case is a uniform flow over an oscillating circular bluff body that can move only in the direction normal to the inlet velocity. The inlet flow velocity is varied systematically such that 90 ≤ Re ≤ 130. Similar structure parameters were used in the numerical simulations of Rosis et al. (2013) and He et al. (2014). The spring is linear, with stiffness k = 5.79 N/m and damping factor c = 0.325 g/s, bluff body mass m = 2.979 g, and bluff body diameter D = 0.16 cm. The fluid is water with viscosity μ = 0.01 g/cms and density is ρ = 1.0 g/cm 3 . In all numerical simulations, the fluid is started from a stationary condition. Figure 5 shows the relative maximum amplitude, A (y,max) , as a function of the Reynolds number, and shows a good agreement between our results and the literature data. The graph also shows the analytical vortex shedding frequency for a stationary circular bluff body.

Results and discussion
In terms of the transverse amplitude, the response characteristics for the square and circular bluff bodies show a considerable discrepancy between the upstream and downstream bluff bodies, respectively. Specifically, the transverse amplitude of the downstream bluff body outperforms that of the isolated bluff body, as shown in Figure 6. In addition, for post-processing, instantaneous time histories of both bluff bodies were documented. The time-averaged response characteristics are evaluated using Equation (11) over a period of 25 oscillations. The response amplitudes of the individual square and circular bluff bodies are added to the graph for understanding and comparison. It is noticeable that the response characteristics for both bluff bodies substantially depend on the spacing ratio (L x /D). The bluff body formation is distributed into four diverse regimes, centered on the behavior of the downstream bluff body.
Regime 1 is categorized as the root of the growing drift in the response amplitude of the downstream bluff body caused by the vigorous coupling of the vortices and the forces, which assists the bluff bodies to considerably oscillate. The results indicate that the downstream bluff body exhibits a considerably higher amplitude in this regime than the upstream one, and the maximum transverse amplitude occurs for the tandem oscillating bluff body cases. In this region, both the tandem square and circular bluff bodies keep shadowing this structure condition A (y,max−upstream) < A (y,max−downstream) . Regime 1 is caused by the highly adjacent locations of the two bluff bodies. In addition, they resemble a single elastic body whose tail end can oscillate to different positions compared with its frontal end. Furthermore, researchers (Adeeb et al., 2018b;Sharma & Eswaran, 2004) reported that for tandem stationary bluff bodies, all incoming flows are faced by the upstream bluff body, whereas the downstream bluff body is entirely enclosed by the upstream bluff body wake. In addition, the stagnation point swings along the leading edge of the bluff bodies, which is a consequence of the oscillatory nature of the gap flow between them.
At the start of regime 1, the upstream bluff body oscillates with an identical magnitude to the corresponding single bluff body. The coupling of the downstream bluff body motion and wake of the upstream bluff body generates an extraordinarily high transverse amplitude (maximum for tandem oscillating simulations) for the downstream square and circular bluff bodies. In addition, for the downstream circular bluff body (L x /D = 1.05), the transverse response amplitude is amplified by 1.5 times compared to those of the single bodies. Overall, the largest amplitude is exhibited by both bluff bodies at the termination of this regime, and the values vary for the square and circular geometries. The response amplitude of the downstream square bluff body is nine times higher (at L x /D = 3.0), whereas for the circular geometry, is twice higher (at L x /D = 1.5) than those of the corresponding single bodies. A remarkable observation is that at the end of this regime, the square bluff body oscillates with an almost equal magnitude to a single circular bluff body. The cause of the shrinking of regime 1 for the circular bluff body is that the flow of the square bluff body is detached at its frontal corners, triggering a long wake length downstream. However, for the circular bluff body, the flow is more attached to the body, causing a small wake length (Adeeb & Sohn, 2021;Sahu & Patnaik, 2011).
As the spacing of the square bluff bodies is further increased (3.0 < L x /D ≤ 4.0), the outcomes decrease, which noticeably reduce the transverse amplitude of the tandem bluff bodies. Although the transverse amplitude of the downstream bluff body decreases to a certain extent, it is still six times higher than that of a single square bluff body. In comparison, for the upstream bluff body, the response approaches that of a single bluff body. In contrast to the square bluff bodies, the upstream and downstream circular bluff bodies (1.5 < L x /D ≤ 3.0) exhibit minimum transverse response amplitudes, which are inferior to the conventional single circular bluff body lock-in amplitude, as shown in Figure 6(b). Furthermore, minimum response values are observed at the termination of this regime. The potential explanation is that these are critical spacing ratios for the geometries; therefore, no vortex street is formed in the gap between the bluff bodies. In addition, the shear layer emitted from the upstream bluff body reattaches to the downstream bluff body. It is noteworthy that various studies (Adeeb et al., 2018b;Sumner, 2010) also show that 4.0D is a critical spacing ratio for tandem square bluff bodies and 3.0D for tandem circular bluff bodies. Thus, in the subsequent discussion, this regime is quantified as a critical spacing or regime 2.
When the bluff bodies are placed further apart, they (square & circular) start to act as two-body systems. The influence of the wake flow of the upstream bluff body becomes dominant, and it is sensed by the downstream bluff body. The response amplitude of the upstream bluff body for both geometries continuously remains close to that of the single bluff body. Moreover, for the downstream bluff body, its amplitude presents a rapid variation. The potential explanation for this is that the vortex starts to roll in front of the downstream bluff body and becomes more distinct, assisting the bluff bodies to achieve higher amplitudes for both geometries again. The span of regime 3 for the square bluff bodies is 4.0 < L x /D ≤ 5.0 and for the circular bluff bodies is 3.0 < L x /D ≤ 4.0. It is worth mentioning that the downstream bluff body (for both geometries) again approaches the maximum transverse response amplitude at the termination of this regime.
By further increasing the streamwise displacement among the tandem square bluff bodies (L x /D > 5.0), a strong distinct wake effect starts to diffuse for the downstream bluff body; nonetheless, it can be neglected. This regime is categorized as regime 4, which can be sub-categorized into two parts: periodic and quasiperiodic. With the increase in the streamwise displacement between the bluff bodies, the transverse amplitude of the downstream bluff bodies starts to decrease, and the instantaneous time history demonstrates a periodic behavior. Figure 6(a) shows that at a large spacing ratio (L x /D > 8.0), the downstream bluff body starts to show a higher response amplitude, which is even equivalent to that in regime 1. However, from the instantaneous time histories, a quasi-periodic behavior is observed at these spacing ratios. Therefore, to interpret the result appropriately, a new graph with the average value of the quasiperiodic pattern is plotted. Figure 6(c) shows that for the downstream square bluff body, the transverse response amplitude is first reduced, and subsequently, it (L x /D > 8.0) becomes nearly constant, which is much lower than that in regime 1. Furthermore, the upstream bluff body always oscillates with an identical magnitude to its corresponding single body. Additionally, the upstream square bluff body exhibits similar quasi-periodic phenomena; however, because its quasi-periodic oscillating amplitude is very small, nearly similar results are obtained.
In contrast to the square bluff bodies, quasi-periodic phenomena are detected first for the circular bluff bodies for streamwise displacements of 5.0 < L x /D ≤ 8.0. Subsequently, to understand correctly, a new graph is plotted for the tandem oscillating circular bluff bodies using the average value of the quasi-periodic pattern, as shown in Figure 6(d). Hence, the results show that during periodic and quasi-periodic parts, the downstream and upstream bluff bodies oscillate at all times with identical magnitudes to their corresponding single bodies. Furthermore, among the stationary tandem bluff bodies, after the critical spacing, at all spacing ratios, the square bluff bodies generate more lift than the circular bluff bodies (Adeeb & Sohn, 2021;Miran & Sohn, 2017). In this regime, the tandem circular upstream and downstream bluff bodies oscillate with identical magnitudes to their corresponding single bodies. However, the upstream square bluff body oscillates with an identical magnitude to the corresponding single square body, whereas the downstream one oscillates similar to the single circular bluff body.

Forces coefficients
Forces play essential roles in explaining the response parameters of tandem bluff bodies. To further explore the influence of the streamwise displacement, time-averaged drag coefficientC D and root mean square lift C (L,rms) coefficients are computed for both the square and circular bluff bodies. Figure 7 shows the time-averaged drag coefficient values for both the upstream and downstream bluff bodies for various spacing ratios between them. For a reasonable comparison, the force values of a single static bluff body and a freely oscillating bluff body are also incorporated in the graph. To acquire the structure natural frequencies, cases of tandem stationary bluff bodies are also simulated, and our results are in good agreement with previous study (Haider & Sohn, 2018). For the tandem stationary bluff bodies, before the critical spacing, the downstream bluff body experiences a negative thrust force (C D ), whereas the upstream bluff body displays a steady decrease in the force. Moreover, for short marginal gaps, both the bluff bodies act as prolonged single bodies. As the streamwise displacement among the bluff bodies is increased, the 'drag-inversion' phenomenon disappears  for the downstream bluff body, and the periodic shedding of the upstream bluff body starts. By a discontinuous jump, almost zero thrust force transforms into a large positive drag force for the corresponding downstream bluff body. The force value for the upstream bluff body becomes equivalent to that of a single bluff body. Moreover, the critical spacing ratios differ for the square and circular bluff bodies.
The results are shown in Figure 7(a,b) for both the square and circular upstream bluff bodies; it can be seen that the drag force always remains near to those of the single oscillating bluff bodies of their respective geometries for large spacing ratios. However, for the square bluff bodies, the values endure during the critical spacing ratio, apart from the spacing ratio, whereas when the spacing ratio becomes very small (L x /D = 1.05), the force value reaches that of a single static square bluff body. In contrast, in the course of the critical spacing ratio, for the circular upstream bluff body, first a steady rise is observed for regime 1, following there is a continuing decrease in the termination of the critical spacing ratio for regime 2. Subsequently, its value reaches that of a single oscillating bluff body, for regimes 3 and 4. The drag-inversion phenomenon is detected for tandem stationary bluff bodies; moreover, because of the oscillating motion of the bluff bodies, drag inversion is not observed for the downstream square and circular bluff bodies.
However, at the start of regime 1, the value is near 0 for the square bluff bodies. As expected, an increase is observed for both geometries in regime 1. After regime 1, for the circular downstream bluff body, the drag value approaches that of a single static bluff body. In comparison, for its counterpart (downstream square bluff body), it approaches 70% of a single stationary bluff body value. For regime 2, as the transverse amplitude is reduced for the downstream (square and circular) bluff bodies, a gradual declining inclination is observed. Subsequently, for regime 3, the well-known discontinuous jump occurs in the values after the critical spacing ratio, as shown in Figure 7(a,b). This is the only point where the values of the circular bluff bodies approach that of a single oscillating bluff body. A possible explanation is that the downstream bluff body begins to face periodic shedding of the upstream bluff body. The drag force magnitudes are significantly different at this spacing ratio; for instance, the square bluff bodies exhibit a magnitude close to that of a single static square bluff body, whereas the circular bluff bodies display an approach to single oscillating bluff bodies. In regime 4, the values start to stabilize and reach a constant value. However, in the quasi-periodic regime, both bluff bodies show a high value. It is noteworthy that a constant trend can be attained by taking the average value of the quasi-periodic pattern (a technique similar to that used for transverse amplitude). For all the spacing ratios of the downstream square and circular bluff bodies, an almost constant value is achieved for the time-averaged periodic and quasi-periodic patterns. Figure 8 shows the results of the root mean square lift, C (L,rms) , for the upstream and downstream bodies for different spacing ratios between the tandem oscillating square and circular bluff bodies. Here significantly different outcomes are attained by the individual bluff bodies. The corresponding lift value intended for a single oscillating square body is less than that of the stationary single square body, whereas a contrast case is observed for the circular bluff bodies. The maximum value is observed for the downstream bluff bodies. The upstream square and circular bluff bodies exhibit lift magnitudes near the respective single oscillating bodies. Moreover, the circular bluff body lift magnitude approaches that of a single stationary circular body: specifically, a small transverse amplitude is observed, which subsequently, gradually recovers. A similar trend is observed for the downstream square and circular bluff bodies owing to their larger magnitude (FIV + WIV phenomena) than the respective upstream bodies. Furthermore, similar maximum values exist for the quasi-periodic regime although if its average value is applied the resultant magnitude approaches the periodic regime. Figures 9 and 10 show the results of the time-averaged velocity deficit of oscillating square and circular bluff bodies, respectively. Overall, five transverse planes are carefully chosen for different arrangements behind the upstream body. Three planes (L g = L x − D; L g = 1/4, 2/4, & 3/4) lie between the upstream and downstream bluff bodies, and two additional planes exist after the downstream body at 1D and 5D streamwise displacements. Moreover, in the first three planes, only the upstream bluff body influences the results, whereas, at the 1D and 5D streamwise displacements, both upstream and downstream bluff bodies influence the result. As described by Assi et al. (2010), in a tandem arrangement, a velocity deficit is sensed by the downstream bluff body compared to the free stream velocity as the time-averaged flow approaches it. In regime 1, as the streamwise displacement is very small, L x /D = 1.5, the planes between the upstream and downstream bluff bodies show almost similar velocity deficits. Furthermore, even after the downstream bluff body (L g = 1D & 5D), an unchanged velocity deficit is observed. However, in contrast with the square bluff body, after the wake region of downstream bluff body velocity deficit effects become to vanish, it progressively attempts to reach the free-stream velocity. In terms of the magnitude of the velocity, the square bluff bodies exhibit greater values as it has a large frontal area that faces the incoming fluid and detaches from its frontal corner (a sharp change in the magnitude).

Velocity deficit
When the bluff bodies are positioned further apart at L x /D = 3.0, both the square and circular geometries display nearly comparable magnitudes of the velocity deficit. For a square bluff body, the velocity deficit area is slightly wider than that of a circular bluff body at all times. The possible cause is that the square and circular upstream and downstream bluff bodies oscillate with similar amplitudes. In both simulations, the velocity deficit is reduced as it moves from the upstream bluff body. Additionally, as the incoming flow travels away from the downstream bluff body (L g = 5D), it starts to recover from the velocity deficit and attempts to approach the free stream velocity.
At this minute as circular bluff bodies are repositioned further apart (L x /D = 4.0), the downstream bluff body comes out of the wake of the upstream bluff body. Consequently, for all five planes, a comparable velocity deficit can be seen (Figure 10(c)). However, the square downstream bluff bodies are still in the wake of the upstream bluff body; therefore, the outcome is similar to that at spacing ratio L x /D = 3.0. Three more graphs are plotted at streamwise displacements L x /D = 6.0, 8.0, and 12.0 to relate the diverse arrangements. For both the square and circular geometries, a large velocity deficit is observed at L g = 5D of the downstream bluff body, and it is noticeable from Figures 16 and 17 that the upstream and downstream vortex energies are combined and produce a large velocity deficit at this position. At the L g = 1D position, the second-largest velocity deficit is observed, and it can be inferred that the downstream bluff body just starts to influence the upstream wake and does not reach the maximum value. Subsequently, in terms of the magnitude, the L g = 1/4 position appears immediately after the upstream bluff body, following which, the free stream velocity decreases the upstream wake effect.
When the bluff bodies are located at L x /D = 12.0, the maximum magnitude is observed for L g = 1D after the downstream bluff body, owing to the combination of the upstream and downstream bluff body wakes. The velocity deficit results start to vanish as the fluid flows in the streamwise direction. Owing to the shape modification, each phenomenon first occurs for a circular bluff body. Consequently, similar phenomena occur at L x /D = 8.0, and circular bluff bodies trail the streamwise spacing ratio L x /D = 12.0 tendency, whereas for the square bluff bodies, this occurs after streamwise spacing ratio L x /D = 8.0.

Spectral analysis
Fourier analysis of the lift force is used in numerous studies to compute the frequency spectrum of a given signal. Figures 11 and 12 present the power spectra of the fluctuating lift force as a function of the frequency ratio for the tandem oscillating square and circular bluff bodies, respectively. The results show that in all simulations of the tandem square and circular bluff bodies, one single frequency peak is always dominant near F r = 1.0 for the individual upstream and downstream bluff bodies. However, even when a second peak arises, the frequency ratio near 1.0 becomes dominant, except for the case of the downstream square body, at which a quasi-periodic amplitude shape is observed. It is worth mentioning that for all spacing ratios of the tandem square bluff bodies, the upstream bluff body always exhibits a smaller magnitude than the downstream bluff body because it has a smaller response amplitude, which is confirmed by fast Fourier analysis. In contrast, for the tandem circular bluff bodies, the downstream circular bluff body continuously persists to dominate until the transverse response amplitude starts to display quasi-periodic WIV   oscillations (L x /D ≤ 4.0). As the streamwise displacement is amplified, the downstream bluff body moves to the quasi-periodic WIV region. In addition, both the upstream and downstream circular bluff bodies start to show the same levels of energy and transverse amplitudes ( Figure 12, L x /D = 7.0, 9.0 and 11.0). For the quasi-periodic region of the tandem oscillating square bluff bodies, the upstream bluff body always oscillates with a very minimal transverse amplitude and without an adequate amount of energy, as supported by Figure  11 (L x /D = 9.0, 11.0). In some cases, multiple dominant frequencies are observed for both the square and circular bluff bodies: one is located adjacent to the natural frequency and the other is situated almost three times to the first frequency. For all simulated cases, an interesting observation is that as the shedding frequency precisely matches the natural frequency, both the square and circular bluff bodies attain high amplitudes.
In Figure 13(a,b), the frequency ratios for the tandem oscillating square and circular bluff bodies are plotted for the entire spacing ratios. For a better understanding, the results of the frequency ratios of the corresponding stationary bluff bodies are also included. A typical synchronization frequency range of (0.8 < F r < 1.2) has been observed in several studies (Goyder, 2002;Haider & Sohn, 2018), and our outcomes are consistent with their results. The results indicate that as the bluff bodies are placed very close to each other, they oscillate close to frequency ratio F r = 1.0. However, as the spacing ratio is increased and moves toward the end region of the critical spacing, they start to show a reduction in the frequency ratio. With further increase in the spacing ratio, the WIV effect starts to dominate for both the square and circular bluff bodies, and they start to oscillate with almost F r = 1.0. One interesting observation for the tandem square bluff bodies is that the quasi-periodic oscillations occurring at the upstream bluff body frequency value remain close to F r = 1.0. In comparison, for the downstream bluff body, the oscillations slightly decrease again, but persist in the typical synchronization frequency range. The potential explanation is that as the downstream square bluff body attains a higher amplitude, its frequency ratio value drops from 1.0, whereas for the circular bluff body, the frequency ratio is close to 1.0.

Phase portrait
To further investigate the dynamic response of the oscillating bluff bodies, a phase portrait graph is plotted between the instantaneous transverse amplitude and the lift force. The results for the square and circular bluff bodies are shown in Figures 14 and 15, respectively. A phase portrait graph provides a significant description of how energy is transferred in fluid and bluff body oscillations. Moreover, the lift force plays a significant role in determining the phase portrait graph profile because it is directly proportional to the pressure distribution. A key observation is that all graphs are symmetric about the origin for both the square and circular bluff bodies at all spacing ratios. This is because both the bluff body shapes are symmetrical, either isolated or in the tandem arrangement. Therefore, the pressure oscillation at the bottom or top surface is also always symmetrical.
The phase portrait graphs also determine whether both the bluff body transverse motion and lift force are in phase or out of phase. The transverse motion and lift force are in phase if the graph is plotted in the first and third quadrants, whereas they are out of phase if graphs are plotted in the second and fourth quadrants. A 45 • phase difference is observed if the graph is plotted inclined to the axes and lies in the first and third quadrants. In contrast, a 135 • phase difference is observed if it is plotted inclined to the axes in the second and fourth quadrants. In a case where it is drawn in all four quadrants, the graph can be divided into four parts. Specifically, first, the bluff body transverse motion and lift force are in phase, which subsequently becomes out of phase; in the third quadrant, they again become in-phase and finally switch to again being out of phase. The converse of these parts can also occur.
For streamwise displacements, L x /D ≤ 8.0, 'except 6D' square bluff bodies, the phase portrait graphs are plotted in the second and fourth quadrants (out of phase).
For L x /D = 6.0 and the reaming streamwise displacement, the phase portrait graphs are plotted in all quadrants, whose significance is already explained above. In contrast to the circular bluff bodies, for all streamwise displacements L x /D > 3.0 and 1.05D, the phase portrait graphs are plotted in the first and third quadrant (in-phase). For a small-time, only the downstream bluff body moves to the opposite quadrants. However, this can be ignored for such a short time. For the reaming  streamwise displacements (1.5D and 2D), graphs are plotted in all four quadrants.
In almost all cases of tandem oscillations, phase portraits of the upstream and downstream bluff bodies are similar, although their amplitude is drastically different. A periodic lock-in condition is observed by the bluff bodies if a single outline cycle is generated, whereas, for the lock-out condition, multiple cycles are noted. The periodic lock-in condition is observed until 5D for square tandem oscillating bluff bodies, where it occurs until 4D for the circular tandem oscillating bluff bodies. With the increase in the spacing ratio for the circular bluff bodies, during the quasi-periodic oscillations, the phase portrait graph shows a very large and overlapping profile compared to that of the square bluff bodies. The profile overlap phenomenon only occurs for the downstream bluff body owing to the WIVs from the upstream bluff body. Similar trends are observed for the square bluff body because of the quasi-periodic oscillations at L x /D > 8.0 (Figure 14 (L x /D = 9.0, 10.0, and 12.0).

Instantaneous vorticity
In line with stationary bluff bodies, time-averaged vortex contours can be plotted and analyzed to gain a better understanding of the wake and the wake length. However, in the case of tandem FIVs and WIVs, both the bluff bodies oscillate; thus, time-averaged vortex contours do not provide any physical significance. The oscillations of the bluff body exchange the fluid portion, which has a velocity deficiency downstream. Consequently, the instantaneous vortex contours are plotted and analyzed. Figures 16 and 17 show the vortex flow patterns for all spacing ratios of the square and circular bluff bodies, respectively. As discussed earlier regarding tandem stationary bluff bodies, owing to critical spacing phenomena, the downstream bluff body is completely covered by the upstream bluff body wake. Similar is observed for the tandem oscillating bluff bodies, and no change occurs in the critical spacing ratio values for the square (4.0D) and circular (3.0D) bluff bodies at the same Reynolds number.
In the range of critical spacing, 1.05 ≤ L x /D ≤ 4.0, the tandem square upstream and downstream bluff bodies shed 2S (single vortex in half oscillation) vortex mode. With a further increase in the spacing ratio, the WIV effects start to appear downstream. A very strong and smooth 2S vortex mode is shed by the upstream square bluff body for the reaming spacing ratios; however, it becomes slightly complex for the downstream bluff body. At L x /D = 5.0, the downstream body sheds a strong 2P (two pairs of vortices in one oscillation) vortex type triggered by the amplified transverse motion at this frequency ratio. For the remaining spacing ratios, the downstream bluff body sheds a P + S(a pair of vortices + a single vortex in one oscillation) vortex mode. Furthermore, for spacing ratios 6.0 ≤ L x /D ≤ 8.0, there is an abrupt change in the vortex profile for the downstream bluff body, causing a long-stretched vortex wake behind the downstream body in this region. For the reaming spacing ratios, however, the shape again alters to a regular vortex shape but with the P + S mode.
In contrast, for the circular bluff bodies, over the critical spacing ratios, both bodies mutually shed a 2S vortex pattern. An exception is spacing ratio L x /D = 1.5, where both bluff bodies achieve the maximum amplitude and the vortex pattern converts to 2P. Similar, to the square bluff bodies, as the circular bluff body leaves the critical spacing ratio, vortex shedding entirely develops above the downstream bluff body. For the reaming spacing ratios, the upstream and downstream circular bluff bodies shed the same vortex mode as the upstream and downstream square bluff bodies (2S for upstream andP + S for downstream), except at L x /D = 4.0 and 5.0. At L x /D = 5.0, a long stretch vortex is generated by the downstream circular bluff body. Additionally, at L x /D = 4.0, the downstream circular body for a second time achieves the maximum transverse amplitude, causing its vortex pattern being converted to the 2P mode.

Conclusions
This study emphasizes the development of an effective numerical solver that simulates FSI phenomena using a computational accelerator. Owing to its parallel scalability, the lattice Boltzmann equation (LBM) is combined with the IBM (to mimic the no-slip boundary condition) as well as with equations of the structural solver of rigid body motion on a CUDA platform. For a Reynolds number of 100 and a frequency ratio of 1.0, various tandem arrangements are investigated for naturally oscillating square and circular bluff bodies. The performance of the GPU accelerator depends on the consumption of the computational resources. The results indicate that the executed GPU code for lid-driven cavity flows updates approximately 500 million lattice sites per second, which is 609 times faster than a CPU code. Furthermore, two more cases are tested: flow over a bluff body (LBM with bounce back) and FSI (IB-LBM) solver, which are 145 and approximately 90 times faster than a serial code, respectively, using a single GPU card. Bottleneck and branch divergence is the foremost drawbacks of the aforementioned; however, it is impossible to avoid them.
Four different flow regimes of bluff body arrangements are recognized based on the response amplitude. (i) A growing trend is detected for the maximum transverse amplitude (L x /D ≤ 3.0 (square) & L x /D ≤ 1.5 (circular)). (ii) A declining development is observed for the transverse amplitude (L x /D ≤ 4.0 (square) & L x /D ≤ 3.0 (circular)). (iii) Increasing inclinations are again observed for the transverse amplitude, with almost identical values of the maximum transverse amplitude at the end of the regime of the corresponding bluff body (L x /D ≤ 5.0 (square) & L x /D ≤ 4.0 (circular)). (iv) Interference of the upstream body wake with the downstream body is observed for all reaming spacing ratios: periodic and quasi-periodic phenomena are observed for both the square and circular bluff bodies at different spacings (L x /D ≤ 12.0 (square and circular). In all simulations, the downstream square and circular bluff bodies always exhibit maximum transverse response amplitudes compared with the respective upstream bodies. The upstream circular bluff body always outperforms the square bluff body for all spacing ratios. A similar trend is observed for the downstream bluff bodies until the spacing ratio of (L x /D ≤ 8.0); however, subsequently, the square bluff bodies enter a quasi-periodic flow where their amplitudes are very comparable to each other.
Similar values for the critical spacing ratio are shown by the static and oscillating tandem square (L x /D ≤ 4.0) and circular bluff bodies (L x /D ≤ 3.0), respectively. In terms of the drag force, for the entire spacing range, the square and circular upstream bodies are adjacent to the respective single oscillating bluff bodies, and the downstream bluff bodies remain near to the respective single static bodies. In the critical spacing ratio, the upstream and downstream circular bodies exhibit a well-known inclining trend even though they oscillate close to the natural frequency. However, the inclining trend is negligible for the square bluff bodies because of the large magnitude of the transverse amplitude. Far from the critical spacing ratios, the lift forces for the upstream and downstream square and circular bluff bodies remain close to the respective single oscillating bluff bodies.
To analyze the time-averaged quantities, velocity deficit graphs are plotted for a comparative study of the square and circular tandem bluff bodies at various spacing ratios. When the streamwise displacement is very small, similar profiles of the velocity deficit are obtained between the gaps. However, at large spacing ratios, the vortex energies of both bodies combine and produce a large velocity deficit behind the downstream square and circular bodies, respectively. Phase portrait graphs are also plotted, which determine whether the bluff body transverse motion and the lift force are in phase or out of phase.
To analyze the instantaneous effect on flow physics, a deep analysis of the vortex structure is conducted in this research. After the critical spacing ratio for the upstream square and circular bluff bodies, the 2S vortex is continually shed in addition to the downstream square and circular bluff body shedding the P + S vortex. An exception is a downstream square body (2P) at L x /D = 5.0 because it attains a nearly maximum amplitude. Moreover, as the circular bluff bodies attain the maximum transverse amplitude (L x /D = 1.5 and 4.0), the downstream bodies start shedding the 2P vortex. In contrast, when the square bluff bodies gain the maximum amplitude, a 2S vortex pattern is shed by the upstream body at critical spacing ratio L x /D = 3.0 and a 2P vortex is shed (downstream body) after critical spacing ratio L x /D = 5.0. For the reaming spacing ratio in the critical range, the tandem square and circular bluff bodies collectively shed a 2S vortex mode. Based on the present results, it is recommended that for tandem configurations, bluff bodies may be arranged such that L x /D = 3.0 and 5.0 (square) & 1.5 and 4.0 (circular) to achieve large force excitations and the resulting transverse amplitudes near the resonance frequency. However, to prevent this, bluff bodies should be place very close to each other with L x /D = 1.05 (square) & 3.0 (circular).

Disclosure statement
No potential conflict of interest was reported by the author(s).