Vibration Performance of A Flow Energy Converter behind Two Side-by-Side Cylinders

Flow-induced vibrations of a ﬂexible cantilever plate, placed in various positions behind two side-by-side cylinders, were computationally investigated to determine optimal location for wake-excited energy harvesters. In the present study, the cylinders of equal diameter D are ﬁxed at center-to-center gap ratio of T / D = 1.7 and immersed in sub-critical ﬂow of Reynold number Re D = 10000. A three-dimensional Navier-Stokes ﬂow solver in an Arbitrary Lagrangian-Eulerian (ALE) description was closely coupled to a non-linear ﬁnite element structural solver that was used to model the dynamics of a composite piezoelectric plate. The cantilever plate was ﬁxed at several positions between 0.5 < x / D < 1.5 and − 0.85 < y / D < 0.85 measured from the center gap between cylinders, and their ﬂow-induced oscillations were compiled and analyzed. Results indicate that ﬂexible plates located at the centerline between the cylinder pairs experience lowest mean in amplitude of oscillation. Maximum overall amplitude in oscillation is predicted when ﬂexible plates are located in the intermediate off-center region downstream of both cylinders. Present ﬁndings indicate potential to further maximize wake-induced energy harvesting plates by exploiting their favourable positioning in the wake region behind two side-by-side cylinders.


Introduction
Flow around single or group of cylinders is a fundamental fluid dynamics problem that have attracted considerable research due to its wide engineering significance. Flow past a cylinder typically involves boundary layer separation/reattachments, free shear layers and vortex shedding, which induces cylinder vibrations and noise generation. These flow dynamics become more complicated when there are interference or interactions with neighbouring cylinders. The understanding behind such flow past multiple cylinders provide important insight as they are found in many branches of engineering applications for example -cooling cores in nuclear reactors, heat exchanger tube bundles, offshore structures, marine risers and pipelines. Among flow past multiple cylinders, the case involving two cylinders has been investigated both numerically and experimentally under several configurations, either in-tandem, side-by-side or staggered arrangements and provide some fundamental understanding of flow or wake interactions [1,2].
Flow characteristics past two side-by-side cylinders are strongly influenced by the gap between cylinders and their diameter-based Reynolds number, Re D . Denoting T as the transverse center-to-center distance between cylinders and D as the cylinder diameter, the flow behaviour around side-by-side cylinders may be classified into three main regimes [3,4] : (i) When T/D < 1.1 − 1.2, the gap between cylinders is small and the flow behaves like passing a single bluff body where vortices are shed alternately from the outer surfaces of the top and bottom cylinders, forming an asymmetric single vortex street. In this singlebluff body regime, the Strouhal number (St = f (2D)/U ∞ ) is approximately 0.2 -similar to a single-bluff body case, but with characteristic length of 2D [5]. (ii) When 1.2 < T/D < 2.0 − 2.2, this intermediate spacing between cylinders is within a critical range where the flow between the gap is bistable and switches direction (flip-flopping) towards either cylinders at irregular intervals (with periods several orders of magnitude larger than the vortex shedding period). As a result, in this biased flow regime, a narrower and wider wake region forms behind either cylinders, corresponding to higher (St n ≈ 0.2 − 0.4) and lower (St w ≈ 0.1 − 0.2) vortex shedding frequencies respectively [6]. Within this regime, the Reynolds number Re may also influence the transition between single and twin vortex streets [7,8]. In their 3-D numerical study, Liu et al. [9] also showed the influence of cylinder inclination (or conversely, flow angle) on formation of these vortex streets. (iii) When T/D > 2.2, the gap is sufficiently large that both cylinders behaves more like isolated bluff bodies resulting in a predominantly symmetric (or anti-phase) parallel wake patterns, that are coupled with a single Strouhal number St of approximately 0.21 [6]. Intermittent in-phase wake patterns may also take place in this symmetric flow regime, but would then synchronizes back to the more stable anti-phase or symmetric pattern. At much higher gaps (T/D > 4.5), any interference associated with proximity of cylinders are negligible and each cylinder behaves as independent bluff bodies with uncoupled flow patterns [3].
Previous studies on flow around bluff bodies with adjacent plates may also provide useful context for the present work. Flow past a single cylinder with plates have received much attention in the past decades, especially for passive flow/vortex control (and subsequently, flow-induced vibration control) of a cylinder (see for example, [10,11,12,13]), understanding of plate dynamics (for instance, [14,15]) and potential for energy harvesting (see for example, [16,17]). Less well explored are flow with plates behind two or multiple cylinders. Furquan et al. [18] computationally investigated vibration response of two flexible splitter plates, each placed behind two side-by-side square cylinders (that are separated by center-to-center distance equal to twice the square edge length L) at Re = 100 under varying reduced velocities (U * = U ∞ / f n L where f n is plate natural frequency). Their results indicate initial anti-phase plate vibration as the plates are attracted towards each other with the accelerating gap flow, before finally synchronizing into an in-phase plate vibration pattern that also exhibits the 'lock-in' phenomena at certain range of U * , as either plates undergo large vibration amplitudes when their response frequency approaches the plate first mode natural frequency. Effect of short splitter plates behind two side-by-side square cylinders (with center-to-center distance equal 3.6 times the square edge length) on their vortex shedding and subsequent sound generation were examined experimentally for Re between 10000-33000 by Octavianty and Asai [19]. In this symmetric flow regime and high Re, generated vortices past the side-by-side square cylinders were highly synchronized and coherent spanwise, allowing effective sound reduction even with shorter splitter plates -in contrast to a single square cylinder where vortices were more three-dimensional and effect of splitter plate was limited. Oruc et al. [20] experimented with splitter plates centrally placed between two side-by-side circular cylinders and examined their effect on the flip-flopping gap flow in the biased flow regime. It was found that at sufficient plate lengths, the asymmetric and bistable wake flow were suppressed, resulting in two symmetric and stable wake patterns.
Demand for clean and sustainable energy that requires minimal human intervention or maintenance (for example, to power wireless sensors in remote areas, electronic devices for structural health monitoring of airborne vehicles or deep ocean structures) have prompted increasing development of various energy harvesting systems. One of a number of prototypes include exploiting wakes behind a single bluff body (for example, [21,17,16,22]) and wakeinduced vibration of cylinders (for example, [23]) to induce oscillations of flexible membranes or plates consisting of piezoelectric materials -converting the fluctuating mechanical strains or energy into electrical energy and have been shown to reach outputs up to 30V [24,25]. Although multiple bluff bodies are ubiquitous in many engineering applications, limited investigation on their potential to excite piezoelectric energy harvesters may be found.
Therefore, in the present study, we aim to extend the current results in the literature by exploring wake-excited piezoelectric energy harvesters placed behind two side-by-side cylinders and also analyze their vibrations at various placements in the wake region with the view of maximizing energy harvesting performance. To that end, we performed computational investigation on unsteady flow past two side-by-side circular cylinders interacting with a thin flexible plate placed at a number of streamwise and crossflow positions behind both cylinders. As the side-by-side cylinders behave like a single bluff body, either in unison or separated when the gap is small or large enough respectively, we considered a center-to-center cylinder spacing of T/D = 1.7 (in the bias flow regime) for the present study, to include effects of multiple bluff body flow interaction or interference. In addition, a moderately high sub-critical flow past cylinders of Re D = 10000 is considered in the present work.

Flow Equations
Although vortex shedding behind a circular cylinder may be initiated at Re D ≈ 49, numerous bluff-body energy harvesting prototypes suggest operations at much higher Reynolds numbers (for example, Re D = 5000 − 40000 for harvesting eel of Allen and Smits [21], Re D = 3000 − 8000 for harvesting concept in Weinstein et al. [26], Re D = 6024 for proposed miniature energy generator in Nguyen et al. [27], Re D in excess of 10000 in the investigation by Yu and Liu [16] and Re D ranging between 3200 to 12000 in the design by Shi et al. [17]). As a result, we consider Re D = 10000 in the present study, which falls in the sub-critical flow regime of flow past circular cylinders. Within this regime, flow is three-dimensional and boundary layers remain laminar prior to separation, which then becomes turbulent in the wake region. In order to better simulate flow separation and turbulence but with more practical computational cost, a recently developed Scale-Resolving or Scale-Adaptive Simulation (SAS) modification on the k − ω Shear Stress Transport (SST) turbulence model proposed by Menter and Egorov [28] was employed. This involves the incompressible, three-dimensional Unsteady Reynolds-Averaged Navier Stokes (URANS) equation for conservation of momentum and the continuity equations, both of which in this case are described in the Arbitrary Lagrangian-Eulerian (ALE) frame of reference in order to accommodate fluid grid movements or deformation, as the grids follow the motion of the interfacing flexible plate [29,30]: where u i , u j and p are the time-averaged fluid velocities, fluid grid velocities and time-averaged fluid pressure respectively (i, j = 1, 2, 3 represent the 3-D cartesian directions). ρ is the constant fluid density, µ denotes the fluid dynamic viscosity and µ t is the turbulent eddy viscosity. In the k − ω SST model, µ t is estimated from the turbulent kinetic (k) energy and specific rate of dissipation (ω) equations. SAS treatment of the k − ω shear stress transport turbulence model estimates a von Kármán length-scale and introduces a corresponding source term into the ω equation, allowing the model to switch between RANS and LES-like simulations [28,31]. Thus, more unsteady turbulent structures are able to be resolved as mesh is refined, but at lower computational costs compared to LES or DNS simulations. Fluid grid velocities are solved using a Laplacian Diffusion model, which smoothly distributes the grid velocities in the whole fluid domain, from a value matching the plate motion near the fluid-plate interface, to zero at the fluid domain boundaries.

Structural Equations
The flexible plate dynamics is described by the unsteady equation of motion in Eq. (3).
where σ ij is the stress tensor and the first term on the right-hand side accounts for the stiffness of the structure, allowing for non-linear large plate deformation and accounting for the effect of induced in-plane tension on bending deflection. In addition, c represents damping ratio, ρ s denotes density of structure material, d represents displacement vector of structure and f represents a time-dependent external force acting on structure (coming from the interfacing fluid pressure and shear). In the present work, we considered neglecting structural damping. Figure 1 shows the computational domain used in the present study, which spans 30D × 16D × 1D in the streamwise (x), crossflow (y) and spanwise (z) directions respectively. The inlet and outlet boundaries are respectively located 10D upstream and 20D downstream of the side-by-side cylinders. As the cylinder centers are spaced T/D = 1.7 apart, we employed a symmetrical height of 8D from either top and bottom boundaries to the middle of the gap between cylinders. In addition, a spanwise width of 1D was considered in the present domain. Studies by Alkishriwi et al. [32] and Shen et al. [33] showed that results with spanwise width of 1D agreed well with experiments. Table 1 summarizes the geometrical parameters employed in the present study for both fluid and structural components.

Computational Model
Boundary conditions prescribed in the computational domain are shown in Figure 1. In addition, both left and right faces of the fluid domain are defined as symmetric conditions and both cylinder surfaces are prescribed as no-slip walls. The fluid faces that interface with the flexible plate are also specified as no-slip walls, implying a dynamic condition that requires matching fluid and plate velocities at these boundaries. For the structural domain, the leading edge of the cantilever flexible plate is fixed, while both left and right side faces of the plate are prescribed as symmetric conditions. Figure 1: Computational domain and boundary conditions. Note that origin of x, y and z coordinate system is located centrally between the cylinder gap.  The physical parameters used for both the fluid and structure are compiled in Tables 2 and 3 respectively. Furthermore, a T/D = 1.7 in the biased flow regime is used. In doing so, we anticipate effects of interference or interactions between flow past side-by-side cylinders may be observed and analyzed. The flexible energy harvesting plate is made of a silicone rubber plate (elastic modulus E = 13.1 MPa, poisson ratio ν = 0.48) with two layers of piezoelectric PVDF sheets (E = 2500 MPa, ν = 0.34) attached along its top and bottom surfaces. The composite properties of this plate are summarized in Table 3. Figure 2 illustrates the unstructured grids employed to discretize the fluid computational domain. A three-dimensional meshing scheme is used where finer meshing is prescribed in the vicinity of the two cylinders and flexible plate. In addition, denser grid points representing the boundary layer meshing is employed close to the surfaces of the cylinders and flexible plate, in order to better resolve their boundary layers. The first grid point normal to both cylinders and flexible plate in the boundary layer mesh were defined 0.001D away from their surfaces, maintaining maximum y+ < 2 in all simulation cases.  A commercial computational fluid dynamics and finite element solver (AN-SYS Workbench 2019 R1) was employed to undertake the flow and structural modelling and coupling. The fluid-structural solver coupling was achieved using a partitioned framework, where the loads (from the fluid solver) were transferred to the structural solver (and applied to the plate) and the resulting plate deformation was feedback to the fluid solver to redefine the fluid domain (and corresponding mesh motion). A number of coupling iterations is typically required within each time step until convergence and a coupling under-relaxation factor of 0.75 was considered for both loads and displacements interchange, to closely couple the fluid-structure iterations. An implicit second order backward euler and Hilbert-Hughes-Taylor (HHT) schemes were used for time-integration of the discretized flow and structural equations respectively. Considering the vortex shedding frequencies of both cylinders and the flip-flopping period of the bistable gap flow that is several orders of magnitude greater than the vortex shedding period, a timestep size of 0.01s was used, which also showed adequate stability in the numerical iterations. The steady-state solution for a given stationary plate condition is used as the initial condition for all the unsteady fluid-structure simulation corresponding to each plate placement cases. Steady state flow condition was typically reached after non-dimensional time tU ∞ /D = 75 (corresponding to at least 15 vortex shedding periods) but were extended to run until tU ∞ /D = 300. Fluid-structure computations were then continued for a further tU ∞ /D = 150. Taking into account initial transients in the fluid-structure solutions, statistical results are presented by time-averaging quantities beyond overall tU ∞ /D = 300.

Simulation Cases
As vortex strength weakens further downstream of a cylinder, we considered a window in the wake region between 0.5 < x/D < 1.5 and −0.85 < y/D < 0.85 to place a flexible plate in the present study. In total, 15 cases were simulated corresponding to locations where leading edge of cantilever plate was fixed, as depicted in Figure 3. These cases are also tabulated in Table 6 highlighting their corresponding streamwise (x) and crossflow (y) coordinate positions.

Validation
A preliminary numerical simulation on two side-by-side cylinders at Re D = 6000 and 10000 was firstly undertaken to verify the present modelling and meshing scheme. Table 4 summarizes comparison of present results against published numerical and experimental data from previous works.   Although C d and C l in the present model are consistent with previous numerical studies (for example, C d = 1.21 and C l = −0.24 for Re D = 6000 at T/D = 1.7 [3] and C d between 1.25 to 1.32 for Re D = 10000 at T/D between 1.5 and 1.75 [34], about which T/D = 1.7 in the present case falls within), the present numerical model overestimated their values in comparison to the experimental results in [37]. However, this difference in C d and C l with experimental results can be seen in all previous numerical studies reported in Table 4. This is potentially due to firstly, the difference in actual experimental to numerical conditions, and secondly, the difficulty in capturing flow separation using turbulence modelling in these sub-critical flow regimes, where the boundary layer remains laminar but the wake region is turbulent. Alternatively, comparison with [38] (albeit for T/D = 1.5) suggests closer agreement with experimental data.
In addition, while the present numerical model show biased gap flow pattern expected at this intermediate gap spacing [2,38], previous numerical studies at Re D = 10000 show symmetric unbiased flow. This may highlight the advantage of the present three-dimensional and Scale-Adaptive Simulation (SAS) turbulence model, which has been shown to better capture the unsteady fluctuations under separated flow conditions than previous URANS turbulence models [31].
Further validation of the present model was performed by undertaking a grid-and time-independence test for flow past two side-by-side cylinders with a stationary plate placed behind the upper cylinder (corresponding to case a1 in Figure 3). A finer grid model (model B), where mesh size was effectively reduced by 1.6 to 2.0 -yielding a total of 226075 nodes in the fluid domain, was constructed. In addition, another simulation using the baseline grid model for case a1 but using a timestep size an order of magnitude lower (ie. model A with ∆t = 0.001s) was undertaken. Table 5 shows small differences between the baseline model results and their corresponding finer grid model and smaller timestep simulation respectively -indicating that results in the present baseline model is sufficiently not sensitive to the size of the mesh and timestep.

Results and Discussion
Overall, 10-20 internal coupling iterations were necessary for fluid-structure coupling convergence at each time step. Computations for the initial stationary plate run may take up to 24 hours, while the following fluid-structure computations may run up to 45 hours, depending on the case simulated. We begin by presenting the oscillation history of the flexible plate when they are placed at each of the 15 locations highlighted in Figure 3. Plate tip displacements (in the y-direction) over time, when plate is positioned along various streamwise locations at each −0.85 < y/D < 0.85, are shown in Figure 4a -4e. Tip oscillations show lowest amplitudes when flexible plate is placed along the centerline (y/D = 0) between the two side-by-side cylinders, compared to other y/D locations. Furthermore, their oscillating amplitude appears to be lowest when plate is closest to the cylinders (x/D = 0.5) before maximizing at x/D = 1.0 and decreases as the plate is placed further downstream, as shown in Figure 4c. In order to better quantify plate vibration behaviour behind the cylinders, the mean y-displacements (y) of the tip was calculated for each case. In addition, as piezoelectric effect is affected by their mechanical strains, which occur irrespective of positive or negative direction in plate deflection, we considered taking their vibration amplitude (A) over time and calculated their root-mean-square (A rms ) as an indication of the overall level in mechanical strains experienced by the piezoelectric beam (and hence, energy harvesting potential), for all 15 cases. Table 6 summarizes these y and A rms values for all the cases simulated in the present study. Inspecting the table, a number of observations may be hy-pothesized in regards to variation in y and A rms with respect to plate placement in the wake region: 1. The mean y-displacements y shows opposing deflections between plate positioned above and below the centerline y/D = 0. At locations immediately behind both cylinders (ie. x/D = 0.5), plate placed above the centerline (y/D > 0) tended to oscillate about a mean position that is deflected downwards (y < 0), while plate placed below the centerline (y/D < 0) tended to oscillate about a mean position that deflects upwards (y>0). However, away from the cylinders (ie. x/D = 1.0.1.5), the opposite occurs, where plate positioned above the centerline appears to oscillate about a mean position that is deflected upwards (ie. y>0) and plate positioned below the centerline appears to deflect more downwards (ie. y<0).  Figure 5 shows the variation in root-mean-square of the oscillation amplitude (A rms ) with respect to plate placement behind two side-by-side cylinders. The contours were generated from the results at the 15 locations in Table 6 using a cubic spline interpolation, giving predicted A rms for the flexible plate if they are located along various (x/D, y/D) locations behind the two cylinders. Thus, Figure 5 indicates potential locations where flexible plate may be placed to undergo high overall vibration amplitudes. This include rear center points on either cylinders or in the intermediate region off-center from the centerline between the cylinders (for placements further downstream). Overall maximum mean amplitude of oscillation is recorded when flexible plate is placed at x/D = 1.5, y/D = −0.425 (ie. case c4), where their mean amplitude is almost twice the mean amplitude on the centerline at similar streamwise (x) position. Figure  6a-6c highlights spanwise vorticity contour plot corresponding to cases with maximum mean amplitude for each streamwise locations x/D = 0.5, 1.0, 1.5. While biased flow regime (with narrow and wide asymmetric wake patterns behind either cylinders) is expected to take place for side-by-side cylinders in the present T/D = 1.7, it is found that this wake pattern still persists in the presence of a flexible plate. Comparing at instances when their respective plate amplitude is maximum, the narrow-wide wake pattern is more pronounced for cases b2 and c2 (ie. Figure 6b and 6c), where plate is placed at y/D = ±0.425, compared to case a5 where plate is placed centrally behind a cylinder at y/D = −0.85. In contrast, the wake pattern for case a3 (corresponding to lowest mean amplitude case located along centerline y/D = 0) in Figure 6d, shows minimal difference in wake size between cylinders, highlighting (i) potential suppression of biased gap flow when flexible plate is centrally positioned and close to both cylinders (similar to previous study in Oruc et al. [20] but with plate length not extending upstream in the present study) and (ii) important role of biased gap flow on plate vibrations behind two side-by-side cylinders. Consequently, wake-excited  Figure 5: Contours representing root-mean-square of plate oscillation amplitudes (A rms ) behind two side-by-side cylinders. energy harvesting plate may be positioned away from the centerline y/D = 0 to maximize their energy output for this type of bluff body configuration.

Conclusion
Wake-excited vibrations of a flexible plate behind two side-by-side cylinders was computationally investigated in their biased flow regime at Re D = 10000. Using a Scale-Adaptive SST turbulence flow model and partitioned coupling to a finite element solver, vibration response of a flexible plate placed at various locations in the wake region were simulated and analyzed. Numerical results indicate placing flexible plate along centerline between both cylinders generated low mean amplitude of oscillations and may potentially suppress biased gap flow if positioned close to both cylinders. It was found that plate vibrations are maximized when placed immediately behind either cylinders or further downstream within some intermediate region off-center from the centerline, where biased asymmetric wake interactions are present. These numerical findings may be exploited to identify favourable placements to maximize wake-excited energy harvesters behind side-by-side circular cylindrical structures. It is remarked that the present study had focused on Re D = 1000, biased flow regime T/D = 1.7 and wake region within 0.5 < x/D < 1.

Funding
This project has received funding from the European Union Horizon 2020 research and innovation programme under the Marie Sklodowska-Curie grant agreement no. 730888.

Abbreviations
The following abbreviations are used in this manuscript: