Numerical Assessment of Flow Energy Harvesting Potential in a Micro-Channel

: A micro-energy harvesting device proposed in the literature was numerically studied. It consists of two bluff bodies in a micro-channel and a flexible diaphragm at its upper wall. Vortex shedding behind bodies induces pressure fluctuation causing vibration of the diaphragm that converts mechanical energy to electrical by means of a piezoelectric membrane. Research on enhancing vortex shedding was justified due to the low power output of the device. The amplitude and frequency of the unsteady pressure fluctuation on the diaphragm were numerically predicted. The vortex shedding severity was mainly assessed in terms of pressure amplitude. The CFD model set-up was described in detail, and appropriate metrics to assess the energy harvesting potential were defined. Several 2D cases were simulated to study the effect of the inlet Reynolds number and channel blockage ratio on the prospective performance of the device. Furthermore, the critical blockage ratio leading to the vortex shedding suppression was sought. A higher inlet velocity for a constant blockage ratio was found to enhance vortex shedding and the pressure drop. Great blockage ratio values but lower than the critical ones seemed to provide great pressure amplitudes at the expense of a moderate pressure drop. There is evidence that the field is fruitful for further research and relevant directions were provided.


Introduction
The development of micro-electronics and wireless communication technology has made wireless sensor networks (WSNs) a very active research field nowadays. WSNs comprise small dimensions sensors and offer important advantages; they can be installed in places where wired connections are not possible, the terrain is inhospitable, or the physical placement of the (self-autonomous and mobile) sensors is challenging. WSNs may prove to be very valuable technology for applications in critical fields, such as national defense, environment monitoring, forest surveillance, healthcare, intelligent buildings, traffic control, industrial process monitoring, target tracking, structural health monitoring, predictive maintenance, etc. [1].
The sensor nodes of WSNs performing data processing, communication, and data transmission are energy-hungry devices. They are equipped with batteries, the energy storage capabilities of which definitely affect the performance of the 'parent' WSN, impose constraints on its sustainable operation, and essentially determine its lifetime. WSNs are usually deployed in places where recharging or replacing their batteries is not feasible. Since the lifetime of a sensor node should usually range from two to ten years, depending on the application, there is no doubt that their batteries, rechargeable or not, are not able to meet this requirement without human presence for replacing or recharging. Although there is ongoing research on increasing the energy density of batteries or reducing the power consumption of WSNs at various levels of their operation (signal processing, operating system, optimization of communication), the most effective method relies on developing techniques that enable the system to repower itself. These are the energy harvesting techniques that generate the required power for each device by harnessing energy from its surroundings. The main idea is that a node could convert energy available in the environment into electrical energy by means of using various conversion schemes and materials. Such an approach eliminates the limitations imposed by the batteries and other finite energy sources of a WSN, so its smooth operation can only be limited due to the failure of any of its components.
The environment has plenty of sources already available to provide unused ambient energy, such as solar, wind, vibration, ocean waves, etc. There is a variety of corresponding methods to harvest this energy. The classification of these methods is mainly based on the different forms of energy involved. The most common methods and sources for harvesting energy are photovoltaic power generation (solar energy), mechanical (piezoelectric, electrostatic, electromagnetic), and dynamic fluid energy harvesting (micro-flows, acoustic, magnetic, hybrid power source). The interested reader can find further information on these topics in relevant review texts [1,2]. However, in the case of powering WSNs, the main disadvantage is that the production scale of the energy destined for miniaturized systems is very small compared to large-scale applications, so millimeter-scale energy harvesting devices are required. Furthermore, when large-scale applications are considered, the related power stations are fixed at a given location, while in small-scale ones, only portable devices are of interest. A remarkable review of micro-scale energy harvesting devices is given in [3]. An important category and portion among these devices utilizes piezoelectric materials, a fact that explains the ongoing research on such materials. Energy harvesting technologies utilizing the piezoelectric phenomenon are reviewed in [4], mainly from the materials' point of view. Among the various energy sources offered for micro-energy harvesting (i.e., mechanical, flow, thermal, solar, electromagnetic, etc.), the present work refers to a case of converting energy from a flow to electricity. The corresponding techniques exploit flow-induced vibrations of piezoelectric membranes either in the form of vortex-induced vibrations or wake-induced vibrations [5]. A review on the latter, i.e., techniques harvesting energy von Karman vortices in airflow to power autonomous sensors, can be found in [6]. A brief but very good and compact literature review on these topics can also be found in the introduction of [7].
In the particular case in which energy is harvested from internal flow in micro-channels, the miniature pneumatic power systems proposed in the literature usually make use either of micro-turbines [8] or bluff bodies [7,9]. The former pose requirements for the precise fabrication of millimeter-scale turbomachinery components, while the latter offer the advantages of simple design, ease of fabrication, and application. In the case of using bluff bodies, one or more of them are appropriately installed into the flow in order to cause enhanced vortex shedding behind them for a wide range of Reynolds numbers. Vortex shedding induces pressure fluctuations that can be exploited by piezoelectric membranes to generate electrical power. Usually, flexible structures that undergo fluidstructure interaction phenomena are used to utilize the piezoelectric effect and harvest flow energy. Various configurations using flexible membranes in conjunction with bluff bodies, both in external and internal flows, have been proposed and assessed in the literature [10,11]. Often, fabricated prototype configurations are experimentally and/or numerically tested to assess their performance [7,12,13]. Some of these studies aim to find the appropriate bluff body shape that produces significant and persistent vortex shedding; such an objective is the same as that used in designing an effective flowmeter [14]. The exploitation of vortex shedding behind a bluff body in the context of designing vortex flowmeters is not new; however, there is continuous research on this. Relevant studies often refer to the effect of bluff body shape on the performance of vortex flowmeters by means of numerical simulations [15,16]. Using multiple bluff bodies in tandem instead of one consists of an interesting perspective toward enhancing the associated vortex shedding. The addition of a second bluff body in order to enhance the repeatability of the vortex shedding pattern has been investigated and verified both numerically [17] and experimentally [7,18]. Experimental studies on the Strouhal number for flows past dual triangulate bluff bodies by varying some geometric parameters have been presented in [19].
In [7], a miniature pneumatic energy-generating device was proposed and tested. This device utilizes one or two bluff bodies in tandem installed in a micro-channel. A flexible diaphragm with a piezoelectric film on it is located above the bodies at the upper channel wall. Pressure fluctuations and unsteady forces induced on the diaphragm due to vortex shedding cause vibrations to it and convert mechanical energy to electrical. Such a device involves ease of fabrication and installation and facilitates miniaturization and massive production, avoiding the need for micro-assembling processes, and could be used in liquid or gas pipeline systems. Its main drawback is the low power output, so further research is required aiming to enhance its performance. The device proposed in [7] was studied by the first author in [20] in order to assess its prospective operation and performance based on numerical simulations. To this end, a 2D CFD model was appropriately set up and the results were compared with corresponding numerical results from [7].
While the device was experimentally tested in [7] by measuring its actual electric power output, in [20], extended use of the CFD model was made in order to study the driving force for the operation of the device, i.e., the pressure fluctuation induced on the diaphragm due to the flow field in the channel, under appropriate assumptions. Various bluff body shapes and configurations were simulated for a fixed inlet flow Reynolds number and channel blockage ratio. Pressure fluctuation characteristics (amplitude and frequency) acting on the center of the diaphragm were numerically predicted. Vortex shedding severity was quantified and assessed in terms of the maximum unsteady pressure fluctuation amplitude. The conclusions of [20] focused on the most effective configuration with respect to body geometry and diaphragm position for a fixed Reynolds number and blockage ratio. Actually, the most effective configuration was found to be the one proposed in [7]. In addition, a greater value of blockage ratio was simulated, and the fact from the literature [7] concerning the corresponding enhancement of pressure fluctuation amplitude was also verified. The work of [7] was continued by the present authors, and new results were presented in a conference paper [21]. In this, a grid independence study was performed, and the effect of the inlet Reynolds number was assessed for a fixed value of the blockage ratio by increasing the inlet flow velocity. The point along the diaphragm experiencing the maximum pressure fluctuation was found and, at that point, it was shown that both the pressure fluctuation amplitude and frequency linearly increase with velocity while the pressure drop due to the bluff bodies in the channel was also calculated.
The present study further extends and completes that of [21], concerning a thorough parametric investigation of the potential performance of the energy harvesting device. The case of [20] with the original Reynolds number and blockage ratio is considered herein to be the baseline case. For the sake of completeness, the set-up of the 2D CFD model [20], the assumptions made, and the definition of appropriate metrics [21] to assess the energy harvesting potential of the device are summarized herein. Detailed grid independence and time-step independence studies were carried out, so grids of sufficient size and timestep values were appropriately selected for the various simulations. The study performed in [21] on the effect of the inlet Reynolds number on the prospective performance of the device is now extended to various channel blockage ratios apart from the baseline one; to this end, several cases were simulated in which the pressure drop was also recorded. Furthermore, the effect of the blockage ratio increase on vortex shedding was studied in detail, and the critical value of BR (causing vortex shedding suppression) was numerically found. The results are presented and discussed, conclusions are drawn, and future research directions are proposed. The main contribution of this work is that it presents a detailed parametric investigation on the potential operation of the device under consideration and provides practical guidelines, from a designer's point of view, on the selection of blockage ratio to enhance pressure amplitude in the expense of a moderate pressure drop.

Set-Up of the CFD Model
The academic version of the ANSYS commercial CFD software FLUENT [22] was used for the numerical simulations. The necessary geometric modeling and grid generation tasks were accomplished by using the relevant modules available in this software. Figure 1 presents the device under consideration where two triangular bodies are located in a flow channel of very small height. A flexible diaphragm is installed above the bodies and causes vibrations to a piezoelectric film connected to it, converting mechanical energy to electrical under the action of the unsteady flow pressure forces; the latter are induced by the vortex shedding occurring behind the bodies. Contrary to most applications where vortex shedding suppression is sought, in such an application, vortex shedding enhancement is sought for a better performance of the device. The flow domain used for the numerical simulations in [20] is shown in Figure 2. The bluff body shape is an isosceles triangle, the base of which faces the incoming flow from the left. The length of the channel is L = 77.06 D, and its height is H = 3.76 D where the width of the bluff body (length of the triangle base) is D = 4.25 mm. The blockage ratio of the channel is defined by the ratio BR = D/H = 0.27. The aspect ratio of each triangle, i.e., the ratio of height to base is 1.95. The flexible diaphragm is located on the upper wall of the channel. Its origin is exactly above the position of the first triangle base. As in [7,20], the distance from the channel inlet to the origin of the diaphragm is 23.53 D = 100 mm. The length of the diaphragm is 8.94 D = 40 mm, and the distance downstream the diaphragm up to the flow outlet is 44.59 D = 160 mm. The whole upper and lower boundaries of the channel are treated as solid walls. The left boundary is the inlet of the flow domain, and the right one is the outlet.

Assumptions Involved in the Simulations
All of the numerical simulations performed in the present work rely on the following assumptions:

•
The flexible diaphragm is considered to be a rigid wall, fluid-structure interaction phenomena are ignored, displacement of the fluid due to the diaphragm motion is ignored, and feedback effects from the diaphragm to the flow are neglected.

•
The diaphragm has small inertia and oscillates with the frequency of vortex shedding; the piezoelectric film is supposed to be strained laterally following the vibrations of the diaphragm and, according to the piezoelectric phenomenon, produce electrical power.

•
Although the actual geometry is three-dimensional, two-dimensional simulations along the symmetry plane of the channel are performed to model the phenomenon.

•
The greater the calculated vortex shedding intensity, the better the expected performance of the device.

Governing Equations and Numerical Solver
In order to simulate the low-speed flow of the air in the micro-channel, the two-dimensional unsteady, incompressible, Reynolds averaged Navier-Stokes (RANS) equations were numerically solved, i.e., the continuity and momentum ones (mean-flow equations). Turbulence was taken into account by means of the realizable variant of the k-ε two-equation model; this variant ensures that only physically realistic viscous stresses will be allowed to arise during the simulations [23]. The mean flow equations were solved by means of FLUENT software [22] based on the SIMPLE pressure correction scheme in the context of the finite volume method on unstructured grids consisting of triangular elements. Second-order spatial accuracy was used for the convective terms of the mean flow equations, while first-order was used for the turbulence model ones. A transient solution of the governing equations was sought by means of a first-order Euler scheme in physical time with a constant time-step, while 20 sub-iterations were applied to converge the solution between two successive time-steps. Inlet velocity was prescribed at the flow inlet (termed as 'velocity inlet' in the software), while the zero-pressure value was set at the outlet (termed as 'pressure outlet'). The no-slip condition was used for the velocity at the walls (channel and bodies). The wall functions were implemented to model velocity profiles at the walls. In particular, standard wall functions [24] were implemented in which the option 'scalable wall functions' offered by the software was activated; the latter ensures that the wall distance employed in wall functions will always be such that y + ≥ y + lim = 11.126. Thus, erroneous modeling of the laminar and buffer boundary layer regions (occurring in the range y + < y + lim) was avoided by effectively shifting the near-wall mesh point to y + = y + lim irrespective of the level of the actual near-wall grid refinement.
For the sake of completeness, the governing equations are written below (where the Einstein convention for summation in repeated indices has been considered, i,j ∈ {1, 2}): Continuity equation Turbulent viscosity (according to Boussinesq approximation) In the above equations, u = (u1, u2) is the mean velocity vector, ρ is the fluid density, p is the pressure, μ is the laminar viscosity, k is the turbulent kinetic energy, ε is the turbulent dissipation rate, Gk is the generation of turbulent kinetic energy, C2 is a constant, while σk and σε are the turbulent Prandtl numbers for k and ε equations, respectively. An essential difference of the realizable variant compared to the standard k-ε model is that the coefficient cμ appearing in Formula (5) for turbulent viscosity is not a constant; it is computed by a relation involving the strain rate and the mean rate-of-rotation tensors. More details on the turbulence model and the constant values can be found in [23].

Description of Case Studies
In the present work, the effect of two parameters was investigated on the energy harvesting potential of the device under consideration. These are the inlet Reynolds number (Re) and the blockage ratio (BR), i.e., the ratio of body width to channel height. To this end, two different studies were carried out; the first concerns the effect of varying the blockage ratio and Reynolds number on the device performance, while the second refers to finding the critical value of the blockage ratio that causes vortex shedding suppression in the channel. In what follows, the cases that had to be simulated in order to accomplish the two studies are described and organized.

Cases to Study the Effect of Blockage Ratio and Reynolds Number
The bluff body width D is the characteristic length of the flow. The Reynolds number (Re) based on the inlet velocity Vin and body width D is Re = ρ·Vin·D/μ. By writing D in terms of the blockage ratio BR as D = BR.·H (H being the channel width), the above formula becomes Re = ρ·Vin ·BR·H/μ. Air was used as the working fluid in all the simulations with a density of ρ = 1.225 kg/m 3 and a dynamic viscosity coefficient of μ = 1.789 × 10 −5 Pa·s. For the baseline case studied in [8] where Vin = 20.7 m/s and BR = 0.27, the Reynolds number was Re = 6024.
Apart from the baseline value of BR = 0.27, 6 more BR values were simulated, i.e., 7 cases in total. For each of these 7 BR values, 7 different Re numbers were simulated (the corresponding cases are named after C1, …, C7), resulting in a total number of 49 simulations with respect to the variation of both Re and BR. In each of these cases, the value of the inlet velocity was set according to the formula Vin = (μ·Re)/(ρ·BR·H). Table 1 summarizes the values of inlet velocity used for the various values of Re and BR. According to the literature, an increase in the blockage ratio generally enhances vortex shedding; however, in confined geometries, there is a critical value of BR that causes suppression of vortex shedding [7]. In order to find this value, a number of simulations were carried out in which the baseline inlet velocity (Vin = 20.7 m/s) was kept constant, while the value of the blockage ratio was gradually increased, starting from a low value up to finding the required critical value.

Fluid Data for the Simulations
The Reynolds number based on channel width H was greater than that based on body width D, i.e., ReH = ρ·Vin·H/μ > ρ·Vin·D/μ = ReD. Thus, according to Table 1, in all cases, Re > 4000 dictates turbulent flow. In order to set boundary conditions for turbulence, the values of turbulence intensity It and turbulent length scale Lt were prescribed at the inlet. Assuming fully developed flow, this was estimated by Lt = 0.07 D, while turbulence intensity was computed as It = 0.16·(ReD) −1/8 [7,20].
A constant time-step was used to march the flow field in physical time. A characteristic time scale of the problem at hand is the so-called 'convective time' defined by TC = D/Vin. As in [7], a value of 5% of Tc was proposed for the physical time-step in the simulations, i.e., Δt = Δtc = (0.05) ·Tc. For the baseline case where D = 0.00425 m and Vin = 20.7 m/s, this value was Δtbsl = 10 −5 s = 0.01 ms. In each case, starting from an ambient initial velocity field, the solution was marched using the corresponding physical time-step for sufficient time to establish periodicity, and then the flow field was allowed to evolve for at least three periods.

Definition of Metrics to Assess the Energy Harvesting Potential of the Device
The assessment of the flow energy harvesting potential in the present work was made on the basis of some appropriately selected and defined metrics. These are quantities of interest, either being indicative for assessing the potential to cause a greater effect on the piezoelectric membrane or characteristic quantities of the flow field, like, for example, the pressure drop along the channel.
Since the time evolution of the pressure acting on the center of the diaphragm exhibits periodicity, its amplitude, i.e., its maximum fluctuation width, was considered to be the basic criterion to comparatively assess the energy harvesting potential of the flow from such a device [8]. To monitor the pressure on the diaphragm, 21 distinct equidistant positions were defined along it on the upper channel wall, denoted by the points P1, P2, …, P21 ( Figure 3). These points serve to evaluate and compare local pressure signals, as well as estimate corresponding amplitude values. The latter is for searching for the maximum energy harvesting potential of the device in the sense that the best performance will be exhibited if the center of the membrane is located at the position experiencing the maximum amplitude of pressure pi(t), i = 1, 2, …, 21. In light of the above, the selected metrics that were computed in each case after establishing flow periodicity are: • The position (point Pi,max among the 21 points) where the maximum pressure fluctuation amplitude Δpi,max is predicted.

•
The value of the maximum pressure fluctuation amplitude Δpi,max and the corresponding non-dimensional quantity (pressure coefficient) Cpi,max = Δpi,max/0.5ρVin 2 .

•
The frequency and the non-dimensional frequency, i.e., Strouhal number St = fD/Vin of the pressure signal at point Pi,max. • The average pressure drop Δpdrop in the duct, calculated as the difference between inlet and outlet average pressures, as well as the same quantity in non-dimensional form, i.e., Cpdrop = Δpdrop/0.5ρVin 2 .
The grid independence study for the baseline case, as well as the assessment of vortex shedding severity for the various inlet Reynolds numbers and blockage ratios to be presented in the following sections have been based on the above-defined metrics.

Grid Generation
In order to numerically simulate the flow in the channel and around the bluff bodies, the flow domain was discretized by means of an unstructured grid consisting of triangular elements. To this end, the grid generation module of the academic version of FLUENT was used. The grid was made denser near the walls in order to better resolve the boundary layer regions. From the various grid metrics available in FLUENT, the parameter named after the 'mesh size' was mainly used to create meshes of varying density.
A series of 7 grids in total were generated to facilitate the grid independence study. The latter is demonstrated herein for the baseline case (Re = 6024, BR = 0.27). Table 2 summarizes the name, size, and corresponding value of the mesh size parameter for each of the seven grids. Figure 4 presents comparative partial views of the grid region near the bodies for the grids corresponding to the mesh size parameter values of 1.0, 0.5, and 0.3.

Grid Independence Study
For the sake of the grid independence study, the baseline case was simulated in all of the grids described in Table 1 (i.e., G1.0 to G0.2). Figure 5a presents the pressure fluctuation amplitude for points P1-P21 as it was predicted in each of the seven grids G1.0 to G0.2. According to the corresponding plots, the position of maximum pressure amplitude along the diaphragm for the first two grids, namely, G1.0 and G0.7, appears to be at positions P8 and P9, respectively, while in all the rest finer grids, it was predicted at point P5, i.e., at a distance 1.88 D = 8 mm from the beginning of the diaphragm (not at the center of it).  Figure 5b presents the average pressure drop in the channel as it was predicted in grids of various sizes. As it becomes evident from this plot, this quantity essentially stops changing after grid G0.3. Figure 6a presents the local pressure evolution in time at point P5 (where the maximum amplitude of pressure fluctuation is predicted) for grids G1.0, G0.7, G0.5, G0.4, and G0.3. It has to be mentioned that the pressure signal curves, after dropping out their transient part until the establishment of periodicity, have been shifted in time for the sake of comparison. According to this figure, it is clear that the various grids from G1.0 to G0.3 exhibit a different resolution of pressure evolution.  Figure 6b presents again the local pressure evolution in time at P5 but for grids G0.3, G0.25, and G0.2. According to it, differences in the pressure evolution between them are practically negligible, i.e., the local pressure evolution at point P5 exhibits negligible variation for grids having a lower size than that of G0.3.
According to the above analysis, grid G0.3 was considered to be sufficient for obtaining a grid-independent solution in the case under consideration, and the size of 0.3 mm was selected to generate all the grids required for the numerical simulations in the present work.
For the sake of completeness, the necessity of using 'scalable wall functions' was assessed by computing the value of y + in several simulations. Thus, for the baseline case in grids of various sizes, the minimum y + was found to be between 0.3-1.1, while the maximum y + was found to be in the range of 17.3-74.0. For the 7 different Reynolds numbers (C1, …, C7), in the case of BR = 0.27 using grid G0.3, the minimum y + was between 0.3-1.1, and the maximum y + was in the range 17.1-27.2. The aforementioned minimum values justify the use of the selected approach.

Time-Step Selection
According to what was prescribed in Section 2.5, the time-step for each simulation was calculated as Δtc = 0.05(D/Vin). In [20] where all the numerical simulations referred to the baseline case (Re = 6400, BR = 0.27), the value computed by this formula, i.e., Δtbsl = 10 −5 s = 0.01 ms, was used. However, in the present study, the inlet velocity in each case was calculated as Vin = (μ.·Re)/(ρ·D), and the above formula for the time-step is written as Δtc = 0.05(ρ·D 2 )/(μ·Re). Furthermore, by taking into account that D = BR.H, the final formula for the time-step as a function of Re and BR is written as Δtc = 0.05(ρ·BR 2 ·H 2 )/(μ·Re), i.e., Δtc depends on BR and Re, so it is different in each case. In particular, the value of Δtc is found in some cases to be lower than Δtbsl (Δtc < Δtbsl); this happens in cases with BR = 0.24 for Re > 4830 (i.e., C3-C7, 5 cases) and in those with BR = 0.27 for Re > 6024 (i.e., C5-C7, 3 cases). For these cases, appropriate comparisons were made in order to assess the effect of the time-step on the resolution of the flow field so as to decide if the use of Δtbsl in any case is sufficient or the corresponding value of Δtc had to be considered instead. The comparison was made in terms of pressure evolution at point P5 where the maximum local pressure fluctuation is predicted in all cases. (Section 4) Figure 7 demonstrates the comparison of the pressure evolution at P5 for cases C5, C6, and C7 with BR = 0.24 and C7 with BR = 0.27, as they are predicted by using either Δt = Δtbsl = 10 −5 s (blue curves) or Δt = Δtc < Δtbsl (red curves). These comparisons reveal that the pressure evolution predicted by the use of the two time-steps Δtc and Δtbsl is different. Furthermore, for the same BR, it seems that the discrepancy between the two curves rather becomes more significant as Re increases (from case C5 to C7). Obviously, in cases where Δtc < Δtbsl, the lower value (i.e., Δtc) is considered to provide more accurate results. Thus, according to the above analysis, the value of the time-step used at each simulation was set as Δt = min{Δtc, Δtbsl}.

Results and Discussion
The operational principle of the device under consideration relies on the severity of vortex shedding. The results to be presented concern the effect of various Reynolds numbers and blockage ratios on its prospective performance, as well as the estimation of the critical blockage ratio leading to vortex shedding suppression. These investigations are made on the basis of the flow metrics defined in Section 2.6.

Effect of Reynolds Number and Blockage Ratio
As mentioned above, the seven cases C1-C7 corresponding to different inlet Reynolds numbers were simulated in grids with characteristics similar to those of G0.3 for each of the seven different BR values. A great part of the results that are presented and discussed below in this section have been produced and presented in the context of [25] Figure 9 presents the variation of pressure amplitude along the diaphragm versus the distance from the beginning of the diaphragm (i.e., for each of the points P1-P21) for various inlet Reynolds numbers in the baseline geometry. According to this, the position of maximum pressure amplitude is found for all cases C1-C7 to be at point P5, i.e., at a distance 8 mm = 1.88 D from the beginning of the diaphragm (at 20% of the diaphragm length from its origin and at the 76% of the distance L = 10.5 mm between the two bodies). Similar plots like that of Figure 9 have been derived for all the other values of blockage ratio, and according to the corresponding results, the maximum pressure amplitude for all cases and all blockage ratios simulated herein was predicted to be at the same location, i.e., at point P5. The authors in [7] state that 'maximum pressure fluctuation is located near where the minimum mean pressure occurs', as it is shown in their plot in Figure 10a reproduced by [7]. More specifically, it seems that the maximum pressure amplitude occurs a little downstream of the location of the minimum mean pressure. This claim motivated the present research to examine if the aforementioned statement concerning the maximum pressure amplitude and mean pressure is also valid for the set of cases simulated herein. The results of this investigation are shown in Figure 10b where the distribution of mean pressure along the diaphragm is plotted. According to this, the position of the minimum mean pressure was found for all cases C1-C7 to be at point P7, i.e., at a distance of 12 mm = 2.82 D from the beginning of the diaphragm (at 30% of the diaphragm length). This seems to validate the statement that the maximum pressure amplitude occurs near the position of the minimum mean pressure; however, herein, in contrast to the results of [7], the maximum pressure amplitude occurs a little upstream of the location of the minimum mean pressure.  Figure 11 presents the pressure evolution at point P5 after periodicity has been established for the baseline geometry and various Reynolds numbers. The results for all cases C1-C7 are distributed and shown in a set of four different plots, each of them containing one or two curves (in order to clearly distinguish them and compare each other). Furthermore, in these four plots, the same range has been used in the vertical axis in order to facilitate cross-comparison among them. So, cases C1-C2 are shown in Figure 11a, C3-C4 in Figure 11b, C5-C6 in Figure 11c, and C7 in Figure 11d. As it can be observed, both the amplitude and fundamental frequency of the pressure evolution at point P5 increase with an increase of Reynolds number.

Effect of Reynolds Number-Various Blockage Ratios
In this subsection, plots of the performance metrics against Re are presented with BR as a parameter. Figure 12a presents the variation of the maximum pressure amplitude at point P5 with respect to the increase of the inlet velocity for various values of BR. A first notice is that the curves in the plot are shifted to the left as the BR increases, i.e., towards lower values of velocity. This is due to the fact that the same range of Reynolds numbers was simulated for each BR, so in each case, the inlet velocity was computed with Vin = (μ·Re)/(ρ·BR·H) (see Section 2.4.2); according to this equation, for the same Re, Vin decreases as BR increases. Referring to Figure 12a, some further remarks can be stated: • An almost linear increase of pressure amplitude with inlet velocity is noticed for all BR; this increase becomes steeper for greater values of BR.

•
The pressure amplitude increases with the increase of BR for the same inlet velocity. Figure 12b presents the corresponding curves in terms of non-dimensional quantities, namely, the coefficient of the maximum pressure amplitude against the Reynolds number. According to this, Cp,max seems to vary slightly for the same BR irrespective of the value of Re, i.e., the effect of Re on Cp,max is very weak. It could be stated that Cp,max is essentially an increasing function of BR only. Figure 13a shows the variation of the fundamental frequency of the pressure signal at point P5 with respect to the inlet velocity for the various BR values. According to it: • An almost linear variation is predicted for all the values of BR; all curves have about the same inclination.

•
For the same inlet velocity, contrary to pressure amplitude, the frequency decreases with the increase of BR.  Figure 13b presents the corresponding curves in terms of non-dimensional quantities, namely, the Strouhal number against the Reynolds number. According to this plot, the value of St for each BR is almost independent of the Re and increases with the increase in BR, i.e., St presents a similar behavior with Cp,max. Figure 14a presents the variation of the average channel pressure drop with the increase of inlet velocity for various BRs. As expected, the pressure drop increases with the square of velocity for a constant BR. For the same inlet velocity, the pressure drop increases for greater values of BR, which is expected, due to the resistance caused in the flow by the greater channel blockage.    Figure 14b presents the corresponding curves in terms of the average pressure drop coefficient Cp,drop against the Re. According to it, Cp,drop slightly varies with Re, however, with a different value for each BR, not monotonically increasing or decreasing with it.

Effect of Blockage Ratio-Various Reynolds Numbers
For the sake of completeness, in this subsection, the results of the previous subsections referring to non-dimensional performance quantities against Re for various BRs are briefly presented in another form, namely, against BR with Re as a parameter. Figure 15a,b shows the plots of the pressure amplitude coefficient Cp,max and Strouhal number, respectively, against BR for various Re; both quantities exhibit an almost linear increase with BR irrespective of Re. Similarly, Figure 16 presents the pressure drop coefficient Cp,drop against BR for various Re; like Cp,max and Strouhal, Cp,drop is almost independent of the Reynolds and increases about linearly with BR except for the final value (this behavior can also be noticed in Figure 14b).

Investigation of Vortex Shedding Suppression Due to Blockage Ratio Increase
Nguyen et al. [7] numerically investigated the effect of the blockage ratio increase on the performance of the device under consideration. To this end, they predicted the mean value and amplitude of the pressure at the center S of the diaphragm for values of BR  higher than the baseline one (i.e., BR > 0.27) for Vin = 20.7 m/s. Figure 17 reproduces their plot from [7] in which the change in the pressure mean value and amplitude at the center of the diaphragm is plotted versus BR. They concluded that by increasing BR, the pressure amplitude is enhanced, attaining a maximum for BR = 0.33, then gradually decreases and becomes 0 for BR = 0.42 where actual vortex shedding suppression is predicted. No other information is provided in [7], for example, on the corresponding behavior of the pressure drop or vortex shedding frequency. Figure 17. Pressure mean value and amplitude at the center of the diaphragm vs. BR (reproduced by [7]).
The authors of the present work attempted to verify the findings of [7], so the BR = 0.42 case was first simulated. Contrary to what the authors claim in [7] reporting suppression of vortex shedding for BR = 0.42, not only is vortex shedding predicted by the present approach at BR = 0.42, but pressure amplitude at the center S of the diaphragm continues to increase significantly beyond that value (for BR > 0.42). This fact motivated the present authors to conduct a thorough study in order to seek the value of BR for which the present CFD model predicts vortex shedding suppression, named after the 'critical' value of BR in what follows.
To this end, a series of appropriate simulations were carried out in which the inlet velocity was kept constant Vin = 20.7 m/s, and the blockage ratio BR was gradually increased while the body aspect ratio value was also kept constant and equal to 1.95 (like in [7]). Each time, a new grid was generated, the corresponding case was simulated, and the results were post-processed up to numerically predict the suppression of vortex shedding. Figure 18 depicts a focused view of some grids used in this study in the region near the bluff bodies for various BR values. In each of these simulations, the mean value and amplitude of pressure at the diaphragm center, the average channel pressure drop, as well as the vortex shedding frequency and Strouhal number were recorded. After several numerical simulations and related trial-and-error efforts, the critical value of BR value was finally found at about BR = 0.662, i.e., significantly higher than that predicted in [3]!  Figure 19a shows the results from the present investigation in a form similar to that of Figure 17, i.e., a change in the pressure mean value and amplitude at the center of the diaphragm versus BR. According to this figure, pressure amplitude increases with the increase in BR, attaining a maximum of BR = 0.655. Then, an abrupt decrease occurs and suddenly becomes 0 for BR = 0.662, i.e., in a very small region of BR increase. This behavior of attaining a maximum and then decreasing to zero is qualitatively the same as that of Figure 17 predicted in [7]. However, these are predicted now in significantly greater values of BR and in a much more pronounced way; the maximum is greater than double, and the corresponding decrease is much more abrupt. Figure 19b provides a focused view of the same plot in the range of the abrupt amplitude decrease, i.e., in the neighborhood of the critical BR value.    For the sake of completeness, Figure 22 demonstrates the predicted pressure signal at the center S of the diaphragm for various values of BR. In particular, the plots in Figure  22a    Discussion on the Estimation of the Critical BR An explanation of the difference between the present results and those of [7] concerning the estimation of the critical BR value is not obvious at all since a very similar numerical approach has been implemented in both works (same commercial software, similar grids, same turbulence model and time-step selection). It has to be mentioned that, concerning the effect of BR on vortex shedding suppression, no further information is provided in [7] for the simulations or any relevant experimental evidence.
For the sake of investigation, the cases BR = 0.42 and BR = 0.55 containing 1 instead of 2 bluff bodies were also simulated. Figure 24a,b presents the evolution of the pressure at the center of the membrane in these cases along with results from the corresponding cases with 2 bodies. According to Figure 24a, the present method predicts vortex shedding for BR = 0.42 even for the case of 1 body and with a predicted amplitude comparable to that of the corresponding case with 2 bodies. According to Figure 24b, vortex shedding has been suppressed for the case with 1 body for BR = 0.55 (but not for the case with the 2 bodies). In parallel with the present work relying on simulations, a corresponding experimental work was conducted in the Department of Naval Architecture of the University of West Attica where the present authors belong. This concerns the fabrication of the device under consideration and the relevant measurements of the produced voltage by the piezoelectric element [26]; in that work, vortex shedding appears in the form of fluctuations in the measured voltage signal. In an effort to have an initial hint on what is the actual behavior of the flow in a case with BR > 0.42, the present authors addressed their question to the group of experimentalists. They asked them to conduct an indicative experiment with Vin = 20.7 m/s and BR = 0.55 in order to check if their measurements are compatible with the occurrence of vortex shedding in such a case. This value of BR is well beyond BR = 0.42, i.e., in a region where the present numerical results clearly predict vortex shedding, while in [7], suppression of vortex shedding has been noticed from the value of BR = 0.42. The desired experiment was performed twice using two piezoelectric membranes of different type, one at a time; with both membranes, it was verified that the measured voltage signal corresponded to a flow exhibiting vortex shedding [27]. For the sake of demonstration, Figure 25a,b depicts the measured voltage in the case with 2 bluff bodies for BR = 0.27 and BR = 0.55, respectively [27]. As can be seen, both cases exhibit similar behavior, that of a flow with vortex shedding. This fact is considered to be positive information for the results of the present study. (Obviously, the joint work that just started with the group of experimentalists must and will definitely be continued).

Further Discussion on the Results
According to the results presented in Section 4.1, the non-dimensional quantities Cp,max, Strouhal number, and Cp,drop all exhibit the same behavior, i.e., are almost constant irrespective of the Re number for the same value of BR. This constant value is an increasing function of BR and could be estimated by averaging for various Res (for example, from Figures 15a,b and 16 for Cp,max, the Strouhal number, and Cp,drop, respectively). Figure 26a presents the variation of the Re-averaged values of Cp,max, St, and Cp,drop versus BR. It can be clearly seen that the increase in BR increases all these three quantities. However, the rate of increase in Cp,max is greater than that of Cp,drop, which again dictates that an appropriate increase in BR may be an efficient way to enhance the vortex shedding effect on the piezoelectric membrane. Figure 26b shows a plot of the ratio Δpmax/ΔpS, i.e., the pressure fluctuation amplitude at P5 to the pressure fluctuation amplitude at the center S of the diaphragm versus the Reynolds number for various values of BR. This ratio slightly decreases with the Reynolds, and for each Re, it increases with BR. Figure 26c shows a plot of the ratio Δpmax/Δpdrop, i.e., the maximum pressure amplitude at P5 to the average channel pressure drop, versus the Reynolds, for various blockage ratios. This ratio seems also to be independent of the Reynolds number and be an increasing function of BR. Figure 26d presents the Re-averaged value for the two above-presented pressure ratios for each value of BR. The following information could be extracted from this figure: • Δpmax is over the BR range from 1.3 to 2.6 times the value of ΔpS (about double in average), and the same is valid for Δpdrop. This means that if the diaphragm was positioned with its center at S at point P5, the achieved pressure amplitude could be multiplied by the corresponding ratio (greater than 100% increase).

•
The value of Δpdrop is of the order of that of ΔpS and definitely lower than Δpmax.
Based on the above statements and according to Figure 20, the greater the BR, the greater the enhancement of ΔpS for BR values lower than the critical one. In particular, for a 44.4% increase of BR (from the baseline value of 0.27 to 0.39), the increase in ΔpS is about 125%. Furthermore, beyond the value of 0.42 that the authors in [7] claim vortex suppression to take place, the increase in ΔpS is much more significant, i.e., for a 41% increase (from 0.39 to 0.55), the increase in ΔpS becomes about 167%. The value of ΔpS for the baseline case according to Figures 9 and 20 is ΔpS = 252 kPa. By further considering the above claim (Figure 26d) that Δpmax~1.5ΔpS (using a moderate estimation), it is concluded that for BR = 0.55, a value of about (252)(1.67)(1.5) = 631.3 kPa could be attained if the center of the membrane was placed at point P5. However, this case corresponds to Re = 12,473, which is outside the range studied in Figure 13b. Furthermore, experimental validation of these important statements becomes absolutely necessary. Anyway, according to all the above, there is evidence that the relevant research field is fruitful and open to further research.

Conclusions-Future Research
A millimeter-scale flow energy harvesting device proposed in the literature [7] was studied numerically. The device contains two bluff bodies installed in a very small flow channel. It exploits vortex shedding behind them to cause oscillations on a flexible diaphragm above them and convert flow energy to electrical one by means of the piezoelectric phenomenon.
In [20], a CFD model was set up for the numerical simulation of this case, and different body shapes and configurations were simulated for a fixed flow Reynolds number and blockage ratio. The achieved vortex shedding severity was assessed in terms of the predicted unsteady pressure fluctuation on the diaphragm in order to find the most efficient configuration, which involved two bodies of triangular section. The study of this configuration was continued initially in [21] and further herein by scrutinizing the CFD model of [20] and performing extended parametric studies to understand the effect of various parameters on the expected performance of the device. In particular, a detailed grid independence study and investigation on the selection of a sufficient time-step for the resolution of the flow phenomena under consideration were performed. The device performance was then numerically assessed in a range of different inlet Reynolds numbers and blockage ratios. Furthermore, the behavior of vortex shedding with respect to the blockage ratio increase was studied in detail, and the critical value of BR (for which vortex shedding suppression occurs) was sought and found.

Conclusions
The conclusions drawn in the present work can be summarized (from a designer's point of view) as: • The maximum pressure amplitude (Δpmax) in all cases occurs at the same position, located upstream of the center of the diaphragm (at a distance of 8 mm from its beginning). Thus, in order to maximize the effect of vortex shedding on the diaphragm, the center of the latter should be placed upstream at the point where the maximum pressure amplitude is predicted.

•
The maximum pressure amplitude increases almost linearly with the inlet velocity for all the values of the blockage ratio (BR); the greater the BR, the more abrupt the increase. Thus, using a greater inlet velocity and greater blockage ratio, a greater maximum pressure fluctuation amplitude can be achieved.

•
The fundamental frequency of the predicted pressure signal at the point where Δpmax occurs increases almost linearly with inlet velocity for all values of BR; the slope of the linear increase remains almost constant for all BR. This frequency slightly decreases with the increase of BR for the same Reynolds number. Since a high frequency is rather desired, maximizing the pressure amplitude (as proposed above) will also lead to a frequency increase.

•
The channel pressure drop (Δpdrop) increases with the square of inlet velocity for all values of BR. For the same inlet velocity, the pressure drop increases with BR. As expected, an increase in pressure amplitude causes an increase in pressure drop.

•
Contrary to similar previous research in the literature [7] for the baseline case, the critical blockage ratio for which vortex shedding suppression occurs was found in the present study to have a significantly greater value, and this seems to be validated by corresponding experiments [27].

•
As a contribution of this work, from a designer point of view and under the prerequisite that these results would be validated by experiments, a great value of BR but lower than its critical one seems to provide a great value of amplitude in the expense of a moderate pressure drop ( Figure 20).

Future Research
In light of the above, aiming to further develop research on the study and design of flow energy piezoelectric harvesting micro-devices, the following directions are proposed for ongoing research and future work:

•
To study of the effect of the distance between the two bluff bodies on the device performance since there is already experimental evidence [26] that a greater distance between the bodies may lead to a greater pressure amplitude.

•
To study the device performance for a particular membrane and attempt to correlate the maximum pressure fluctuation amplitude predicted by CFD with the measured electric power from the corresponding experiments, i.e., to extract the operational curve of the device.

•
To perform design optimization studies with respect to characteristic geometric quantities (location and distance between the two bodies, location of the diaphragm, channel blockage ratio, etc.) for maximum performance. A stochastic-based approach can be implemented (e.g., genetic algorithms), which would utilize either the present CFD solver or any reduced-order model of the phenomenon. The solution of a multiobjective problem could be sought, e.g., maximization of Δpmax with minimization of Δpdrop. Furthermore, this could be a constrained problem, e.g., by requiring the frequency to be near the resonant frequency of the membrane.

•
To model the phenomenon more accurately, like, for example, to compare 3D against 2D simulations and/or model membrane dynamics and to consider fluid-structure interaction in the simulations.