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

Abstract

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.

1. 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 fluid-structure 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 time-step 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.

2. 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.

2.1. Flow Domain Geometry

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.

2.2. 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.

2.3. 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

(1)$\frac{\partial u_j}{\partial x_j}=0$

Momentum equation

(2)$\rho \frac{\partial u_i}{\partial t}+\rho u_j\frac{\partial u_i}{\partial x_j}=-\frac{\partial p}{\partial x_i}+\left(\mu +{\mu }_t\right)\frac{{\partial }^2{u}_i}{\partial x_j\partial x_j}$

k-equation

(3)$\frac{\partial (\rho k)}{\partial t}+\frac{\partial \left(\rho k{u}_j\right)}{\partial x_j}=\frac{\partial }{\partial x_j}\left[\left(\mu +\frac{{\mu }_t}{{\sigma }_k}\right)\frac{\partial k}{\partial x_j}\right]+G_k-\rho \epsilon$

ε-equation

(4)$$\begin{gathered} \frac{\partial (\rho \epsilon )}{\partial t}+\frac{\partial \left(\rho \epsilon u_j\right)}{\partial x_j}=\frac{\partial }{\partial x_j}\left[\left(\mu +\frac{{\mu }_t}{{\sigma }_k}\right)\frac{\partial \epsilon }{\partial x_j}\right]+\rho C_1\epsilon \sqrt{2S_{ij}{S}_{ij}}-\rho C_2\frac{{\epsilon }^2}{k+\sqrt{\mu \epsilon /\rho }} \\ {C}_1=\max \left\{0.43,\frac{\eta }{\eta +5}\right\},\eta =\frac{k}{\epsilon }\sqrt{2S_{ij}{S}_{ij}},S_{ij}=\frac{1}{2}\left(\frac{\partial u_i}{\partial x_j}+\frac{\partial u_j}{\partial x_i}\right),G_k=2{\mu }_t{S}_{ij}{S}_{ij} \end{gathered}$$

Turbulent viscosity (according to Boussinesq approximation)

(5)${\mu }_t=\rho c_{\mu }\frac{k^2}{\epsilon }$

Standard wall functions $u^{+}=\{\begin{array}{l}\ln \left(9.793y^{+}\right)/0.4187,\hspace{1em}{y}^{+}\ge y_{\lim }^{+}\\ y^{+}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em},\hspace{1em}{y}^{+}<y_{\lim }^{+}\end{array}$.

Scalable wall functions $y^{+}=\max \left\{y^{+},y_{\lim }^{+}\right\}$.

In the above equations, u = (u₁, u₂) 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, G_k is the generation of turbulent kinetic energy, C₂ 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].

2.4. 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.

2.4.1. 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 V_in and body width D is Re = ρ·V_in·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 = ρ·V_in ·BR·H/μ. Air was used as the working fluid in all the simulations with a density of ρ = 1.225 kg/m³ and a dynamic viscosity coefficient of μ = 1.789 × 10⁻⁵ Pa·s. For the baseline case studied in [8] where V_in = 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 V_in = (μ·Re)/(ρ·BR·H). Table 1 summarizes the values of inlet velocity used for the various values of Re and BR.

2.4.2. Cases to Find the Value of Blockage Ratio Causing Suppression of Vortex Shedding

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 (V_in = 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.

2.5. Fluid Data for the Simulations

The Reynolds number based on channel width H was greater than that based on body width D, i.e., Re_H = ρ·V_in·H/μ > ρ·V_in·D/μ = Re_D. 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 I_t and turbulent length scale L_t were prescribed at the inlet. Assuming fully developed flow, this was estimated by L_t = 0.07 D, while turbulence intensity was computed as I_t = 0.16·(Re_D)^−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 T_C = D/V_in. As in [7], a value of 5% of T_c was proposed for the physical time-step in the simulations, i.e., Δt = Δt_c = (0.05) ·T_c. For the baseline case where D = 0.00425 m and V_in = 20.7 m/s, this value was Δt_bsl = 10⁻⁵ 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.

2.6. 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 p_i(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 P_i,max among the 21 points) where the maximum pressure fluctuation amplitude Δp_i,max is predicted.

• The value of the maximum pressure fluctuation amplitude Δp_i,max and the corresponding non-dimensional quantity (pressure coefficient) Cp_i,max = Δp_i,max/0.5ρV_in².

• The frequency and the non-dimensional frequency, i.e., Strouhal number St = fD/V_in of the pressure signal at point P_i,max.

• The average pressure drop Δp_drop 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., Cp_drop = Δp_drop/0.5ρV_in².

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.

3. Selection of Grid Size and Time-Step

3.1. 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.

3.2. 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., G_1.0 to G_0.2).

Figure 5a presents the pressure fluctuation amplitude for points P1–P21 as it was predicted in each of the seven grids G_1.0 to G_0.2. According to the corresponding plots, the position of maximum pressure amplitude along the diaphragm for the first two grids, namely, G_1.0 and G_0.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 G_0.3.

Figure 6a presents the local pressure evolution in time at point P5 (where the maximum amplitude of pressure fluctuation is predicted) for grids G_1.0, G_0.7, G_0.5, G_0.4, and G_0.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 G_1.0 to G_0.3 exhibit a different resolution of pressure evolution.

Figure 6b presents again the local pressure evolution in time at P5 but for grids G_0.3, G0.₂₅, and G_0.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 G_0.3.

According to the above analysis, grid G_0.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 G_0.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.

3.3. Time-Step Selection

According to what was prescribed in Section 2.5, the time-step for each simulation was calculated as Δt_c = 0.05(D/V_in). 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., Δt_bsl = 10⁻⁵ s = 0.01 ms, was used. However, in the present study, the inlet velocity in each case was calculated as V_in = (μ.·Re)/(ρ·D), and the above formula for the time-step is written as Δt_c = 0.05(ρ·D²)/(μ·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 Δt_c = 0.05(ρ·BR²·H²)/(μ·Re), i.e., Δt_c depends on BR and Re, so it is different in each case. In particular, the value of Δt_c is found in some cases to be lower than Δt_bsl (Δt_c < Δt_bsl); 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 Δt_bsl in any case is sufficient or the corresponding value of Δt_c 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 = Δt_bsl = 10⁻⁵ s (blue curves) or Δt = Δt_c < Δt_bsl (red curves). These comparisons reveal that the pressure evolution predicted by the use of the two time-steps Δt_c and Δt_bsl 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 Δt_c < Δt_bsl, the lower value (i.e., Δt_c) 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{Δt_c, Δt_bsl}.

4. 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.

4.1. 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 G_0.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].

4.1.1. Effect of Reynolds Number—Baseline Geometry

Some representative pictures of the flow field are presented first for the sake of demonstration. These refer to the baseline geometry (BR = 0.27). In particular, Figure 8a shows iso-velocity contours predicted at time 8 ms for case C4 (V_in = 20.7 m/s, Re = 6024), while Figure 8b shows corresponding results for case C7 (V_in = 26.7 m/s, Re = 7770). Obviously, due to the greater inlet velocity, greater velocities are attained in the case of C7 above and below the bluff bodies (red regions) where the flow is accelerated due to the narrow flow passage.

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.

4.1.2. 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 V_in = (μ·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, C_p,max seems to vary slightly for the same BR irrespective of the value of Re, i.e., the effect of Re on C_p,max is very weak. It could be stated that C_p,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 C_p,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 C_p,drop against the Re. According to it, C_p,drop slightly varies with Re, however, with a different value for each BR, not monotonically increasing or decreasing with it.

4.1.3. 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 C_p,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 C_p,drop against BR for various Re; like C_p,max and Strouhal, C_p,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).

4.2. 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 V_in = 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.

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 V_in = 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.

Figure 20 presents in the same diagram, the pressure amplitude at S with the increase of BR, along with the corresponding change in the channel average pressure drop. The comparative view of these two curves, namely, that of the pressure amplitude (for which a high value is desired) and that of the pressure drop (that a high value is unwanted) may be particularly useful to the designer of the device. Their comparison shows that, in the case under consideration, the increase in pressure amplitude due to the increase of the body width comes at the cost of about a similar pressure drop up to the value of BR = 0.42. Then, the two curves diverge, and they come close again at BR = 0.65 where both the amplitude and pressure drop are maximized. From a designer’s point of view and provided that these results would be validated by experiments, the value of BR = 0.60 seems to provide a significant increase in amplitude at the expense of a moderate pressure drop.

Figure 21a,b presents the corresponding plots of frequency and the Strouhal number versus BR, respectively (always referring to the pressure signal at the center S of the diaphragm). Both quantities follow the behavior described for the pressure amplitude in Figure 20; the predicted vortex shedding beyond the value of BR = 0.42 is characterized by an increase in the frequency up to BR = 0.660 where it is maximized and then exhibits a steep decrease.

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 refer to the values of BR = 0.65, 0.655, 0.660, 0.661, 0.662, 0.663, while those in Figure 22b refer to BR = 0.661, 0.662, 0.663. As it becomes evident from these plots, the value of BR = 0.662 is the first one for which the pressure essentially stops fluctuating.

Figure 23 demonstrates iso-surfaces of instant velocity at a time 25 ms after the beginning of the flow solution for some BR values. Vortex shedding seems to be suppressed for BR = 0.662.

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 V_in = 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).

4.3. Further Discussion on the Results

According to the results presented in Section 4.1, the non-dimensional quantities C_p,max, Strouhal number, and C_p,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 Figure 15a,b and Figure 16 for C_p,max, the Strouhal number, and C_p,drop, respectively).

Figure 26a presents the variation of the Re-averaged values of C_p,max, St, and C_p,drop versus BR. It can be clearly seen that the increase in BR increases all these three quantities. However, the rate of increase in C_p,max is greater than that of C_p,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 Δp_max/Δp_S, 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 Δp_max/Δp_drop, 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:

• Δp_max is over the BR range from 1.3 to 2.6 times the value of Δp_S (about double in average), and the same is valid for Δp_drop. 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 Δp_drop is of the order of that of Δp_S and definitely lower than Δp_max.

Based on the above statements and according to Figure 20, the greater the BR, the greater the enhancement of Δp_S 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 Δp_S is about 125%. Furthermore, beyond the value of 0.42 that the authors in [7] claim vortex suppression to take place, the increase in Δp_S is much more significant, i.e., for a 41% increase (from 0.39 to 0.55), the increase in Δp_S becomes about 167%. The value of Δp_S for the baseline case according to Figure 9 and Figure 20 is Δp_S = 252 kPa. By further considering the above claim (Figure 26d) that Δp_max~1.5Δp_S (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.

5. 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.

5.1. Conclusions

The conclusions drawn in the present work can be summarized (from a designer’s point of view) as:

• The maximum pressure amplitude (Δp_max) 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 Δp_max 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 (Δp_drop) 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).

5.2. 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 Δp_max with minimization of Δp_drop. 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.

Author Contributions

Conceptualization, methodology, investigation, writing—original draft, writing—review and editing, preparation, supervision; D.G.K., software, validation, investigation, resources, visualization, writing—review and editing; M.V.N.B. All authors have read and agreed to the published version of the manuscript.

Data Availability Statement

Data sharing is not applicable.

Acknowledgments

D.N. Pagonis and graduate student I. Matsoukas, Department of Naval Engineering, School of Engineering, University of West Attica are gratefully acknowledged for conducting the experiment to validate the occurrence of vortex shedding in the case of BR = 0.55.

Conflicts of Interest

The authors declare no conflict of interest.

Footnotes

Disclaimer/Publisher’s Note: The statements, opinions and data contained in all publications are solely those of the individual author(s) and contributor(s) and not of MDPI and/or the editor(s). MDPI and/or the editor(s) disclaim responsibility for any injury to people or property resulting from any ideas, methods, instructions or products referred to in the content.

Figures and Tables

Figure 1. Miniature energy harvesting device using two triangles in tandem (from [3]): (a) configuration, (b) operational principle.

Figure 2. View of the computational flow domain used for the numerical simulations (not to scale).

Figure 3. Positions P1, P2, …, P21 along the membrane.

Figure 4. Grids corresponding to mesh size values of 1 mm (left), 0.5 mm (middle), and 0.3 mm (right).

Figure 5. (a) Position of maximum pressure fluctuation amplitude along diaphragm for various grid sizes. (b) Average pressure drop predicted for various size grids.

Figure 6. Local pressure evolution at point P5 for grids (a) G1.0 to G0.3 and (b) G0.3 to G0.2.

Figure 7. Pressure evolution at point P5 for cases (a) C5 with BR = 0.24, (b) C6 with BR = 0.24, (c) C7 with BR = 0.24, and (d) C7 with BR = 0.27. The blue curves correspond to results obtained with Δt = Δtbsl = 10−5 s, while the red ones to results obtained with Δt = Δtc < Δtbsl.

Figure 8. Iso-velocity contours predicted at 8 ms in the baseline geometry (BR = 0.27) for cases (a) C4 (Re = 6024, Vin = 20.7 m/s) and (b) Re = 7770 (Vin = 26.7 m/s).

Figure 9. Pressure fluctuation amplitude versus distance from the beginning of the diaphragm for baseline geometry (BR = 0.27) and various inlet velocities (Re).

Figure 10. (a) Pressure amplitude and mean pressure versus streamwise coordinate for the baseline case (from [3]). (b) Mean pressure versus distance from the beginning of diaphragm for baseline BR and various Re (present results).

Figure 11. Comparative pressure evolution for the baseline geometry at point P5 for cases (a) C1 (blue)–C2 (orange), (b) C3 (grey)–C4 (yellow), (c) C5 (cyan)–C6 (green), and (d) C7 (the same range has been used in the vertical axis).

Figure 12. (a) Pressure amplitude vs. inlet velocity and (b) pressure amplitude coefficient vs. Reynolds number.

Figure 13. (a) Frequency vs. inlet velocity and (b) Strouhal vs. Reynolds number for various BR.

Figure 14. (a) Average pressure drop vs. inlet velocity and (b) average pressure drop coefficient vs. Re for various values of BR.

Figure 15. (a) Pressure amplitude and (b) Strouhal number against BR for various Re.

Figure 16. Average pressure drop coefficient vs. BR for various Re.

Figure 17. Pressure mean value and amplitude at the center of the diaphragm vs. BR (reproduced by [7]).

Figure 18. Focused view of grids near the bluff bodies for various BR values in increasing order.

Figure 19. (a) Mean value (in blue color) and amplitude (red bars) of pressure at the center of the diaphragm vs. BR by the present method. (b) Focused view in the range near the critical BR value.

Figure 20. Pressure amplitude on the center S of the diaphragm and average channel pressure drop as functions of BR for constant inlet velocity.

Figure 21. (a) Frequency and (b) Strouhal number of the pressure signal at the center of the diaphragm as functions of BR for constant inlet velocity.

Figure 22. Pressure signal at the center S of the diaphragm. (a) For values of BR = 0.65, 0.655, 0.660, 0.661, 0.662, 0.663. (b) For values of BR = 0.661, 0.662, 0.663.

Figure 23. Instant iso-surfaces of velocity 25 ms after beginning the flow solution. Clearly, vortex shedding seems to be suppressed for BR = 0.662.

Figure 24. Pressure evolution at the center of the membrane for the cases with 1 or 2 bluff bodies for (a) BR = 0.42 and (b) BR = 0.55.

Figure 25. Measured voltage in case with 2 bluff bodies for (a) BR = 0.27 and (b) BR = 0.55 [27].

Figure 26. (a) Re-averaged values of Cp,max, Strouhal, and Cp,drop versus BR. (b) Ratio Δpmax/ΔpS versus Re for various BR. (c) Ratio Δpmax/Δpdrop versus Re for various BR. (d) Re-averaged values of Δpmax/ΔpS and Δpmax/Δpdrop versus BR.

Table 1

Inlet velocity used for the various Reynolds numbers and blockage ratios (data for baseline case are in bold).

Case	C1	C2	C3	C4	C5	C6	C7
Re	4278	4860	5442	6024	6606	7188	7770	BR
	16.3	18.5	20.7	22.9	25.1	27.3	29.6	0.24
	14.7	16.7	18.7	20.7	22.7	24.7	26.7	0.27
	13.0	14.8	16.6	18.3	20.1	21.9	23.6	0.30
Vin (m/s)	11.8	13.4	15.1	16.7	18.3	19.9	21.5	0.33
	10.8	12.3	13.8	15.3	16.7	18.2	19.7	0.36
	10.0	11.4	12.7	14.1	15.5	16.8	18.2	0.39
	9.3	10.6	11.8	13.1	14.4	15.6	16.9	0.42

Table 2

Mesh sizes used in the grid independence study.

Name	G1.0	G0.7	G0.5	G0.4	G0.3	G0.25	G0.2
Mesh size (mm)	1.00	0.70	0.50	0.40	0.30	0.25	0.20
# Nodes	6415	12,623	21,433	30,490	56,529	87,216	105,783
# Elements	12,086	24,230	41,491	59,288	110,840	171,737	208,238

References

1. Tan, Y.K.; Panda, S.K. Review of Energy Harvesting Technologies for Sustainable Wireless Sensor Networks. Sustainable Wireless Sensor Networks; 2nd ed. Kheng, T.Y.; Sheah, W. IntechOpen: London, UK, 2010; pp. 15-43. Available online: https://www.intechopen.com/chapters/12418 (accessed on 17 July 2023).

2. Zhou, G.; Huang, L.; Li, W.; Zhu, Z. Harvesting ambient environmental energy for wireless sensor networks: A survey. J. Sens.; 2014; 2014, 815467. [DOI: https://dx.doi.org/10.1155/2014/815467]

3. Selvan, K.V.; Mohamed Ali, M.S. Micro-scale energy harvesting devices: Review of methodological performances in the last decade. Renew. Sustain. Energy Rev.; 2016; 54, pp. 1035-1047. [DOI: https://dx.doi.org/10.1016/j.rser.2015.10.046]

4. Safaei, M.; Sodano, H.A.; Anton, S.R. A review of energy harvesting using piezoelectric materials: State-of-the-art a decade later (2008–2018). Smart Mater. Struct.; 2019; 28, 113001. [DOI: https://dx.doi.org/10.1088/1361-665X/ab36e4]

5. Zhang, M.; Wang, J. Experimental study on piezoelectric energy harvesting from vortex-induced vibrations and wake-induced vibrations. J. Sens.; 2016; 2016, 2673292. [DOI: https://dx.doi.org/10.1155/2016/2673292]

6. Demori, M.; Ferrari, M.; Bonzanini, A.; Poesio, P.; Ferrari, V. Autonomous sensors powered by energy harvesting from von Karman vortices in airflow. Sensors; 2017; 17, 2100. [DOI: https://dx.doi.org/10.3390/s17092100] [PubMed: https://www.ncbi.nlm.nih.gov/pubmed/28902139]

7. Nguyen, H.-D.T.; Pham, H.-T.; Wang, D.-A. A miniature pneumatic energy generator using Karman vortex street. J. Wind Eng. Ind. Aerodyn.; 2013; 116, pp. 40-48. [DOI: https://dx.doi.org/10.1016/j.jweia.2013.03.002]

8. Lyshevski, S.E. High-power density miniscale power generation and energy harvesting systems. Energy Conv. Manag.; 2011; 52, pp. 46-52. [DOI: https://dx.doi.org/10.1016/j.enconman.2010.06.030]

9. Wang, D.-A.; Chiu, C.-Y.; Pham, H.-T. Electromagnetic energy harvesting from vibrations induced by Karman vortex street. Mechatronis; 2012; 22, pp. 746-756. [DOI: https://dx.doi.org/10.1016/j.mechatronics.2012.03.005]

10. Akaydin, H.D.; Elvin, N.; Andreopoulos, Y. Wake of a cylinder: A paradigm for energy harvesting with piezoelectric materials. Exp. Fluids; 2010; 49, pp. 291-304. [DOI: https://dx.doi.org/10.1007/s00348-010-0871-7]

11. Sanchez-Sanz, M.; Fernandez, B.; Velazquez, A. Energy-harvesting microresonator based on the forces generated by the Karman street around a rectangular prism. J. Microelectromech. Syst.; 2009; 18, pp. 449-457. [DOI: https://dx.doi.org/10.1109/JMEMS.2009.2013395]

12. Hu, Y.; Yang, B.; Chen, X.; Wang, X.; Liu, J. Modeling and experimental study of a piezoelectric energy harvester from vortex shedding-induced vibration. Energy Conv. Manag.; 2018; 162, pp. 145-158. [DOI: https://dx.doi.org/10.1016/j.enconman.2018.02.026]

13. Gharghani, F.M.; Bijarchi, M.A.; Mohammadi, O.; Shafii, M.B. An experimental investigation into a novel small-scale device for energy harvesting using vortex-induced vibration. Int. J. Low-Carbon Technol.; 2021; 16, pp. 317-325. [DOI: https://dx.doi.org/10.1093/ijlct/ctaa060]

14. Miau, J.J.; Liu, T.W. Vortex flowmeter designed with wall pressure measurement. Rev. Sci. Instrum.; 1990; 61, pp. 2676-2681. [DOI: https://dx.doi.org/10.1063/1.1141859]

15. Gandhi, B.K.; Singh, S.N.; Seshadri, V.; Singh, J. Effect of bluff body shape on cortex flow meter performance. Int. J. Eng. Mater. Sci.; 2004; 11, pp. 378-384.

16. Rzasa, M.R.; Czapla-Nielacna, B. Analysis of the Influence of the Vortex Shedder Shape on the Metrological Properties of the Vortex Flow Meter. Sensors; 2017; 21, 4697. [DOI: https://dx.doi.org/10.3390/s21144697] [PubMed: https://www.ncbi.nlm.nih.gov/pubmed/34300437]

17. Jan, Y.-J.; Sheu, T.W.-H. A numerical confirmation of the dual body vortex flowmeter. Comput. Fluids; 2004; 33, pp. 1157-1174. [DOI: https://dx.doi.org/10.1016/j.compfluid.2003.09.004]

18. Peng, J.; Fu, X.; Chen, Y. Flow measurement by a new type vortex flowmeter of dual triangulate bluff body. Sens. Actuators A Phys.; 2004; 115, pp. 53-59. [DOI: https://dx.doi.org/10.1016/j.sna.2004.03.020]

19. Peng, J.; Fu, X.; Chen, Y. Experimental investigations of Strouhal number for flows past dual triangulate bluff bodies. Flow Meas. Instrum.; 2008; 19, pp. 350-357. [DOI: https://dx.doi.org/10.1016/j.flowmeasinst.2008.05.002]

20. Koubogiannis, D.G. Parametric CFD study of micro-energy harvesting in a flow channel exploiting vortex shedding. Open Eng.; 2016; 6, pp. 135-144. [DOI: https://dx.doi.org/10.1515/eng-2016-0016]

21. Benetatos, M.V.; Koubogiannis, D.G. Numerical Investigation of the Energy Harvesting Potential from a Micro-channel Flow. Proceedings of the 10th International Conference from “Scientific Computing to Computational Engineering” (10th IC-SCCE); Athens, Greece, 6–9 July 2022.

22. Ansys Fluent, Academic Version. Available online: https://www.ansys.com/academic (accessed on 24 June 2023).

23. Shih, T.-H.; Liou, W.W.; Shabbir, A.; Yang, Z.; Zhu, J. A New k-ε Eddy-Viscosity Model for High Reynolds Number Turbulent Flows—Model Development and Validation. Comput. Fluids; 1995; 24, pp. 227-238. [DOI: https://dx.doi.org/10.1016/0045-7930(94)00032-T]

24. Launder, B.E.; Spalding, D.B. The Numerical Computation of Turbulent Flows. Comput. Methods Appl. Mech. Eng.; 1974; 3, pp. 269-289. [DOI: https://dx.doi.org/10.1016/0045-7825(74)90029-2]

25. Benetatos, M.V. Assessment of Flow Energy Harvesting in a Duct by Means of CFD Numerical Simulations. Diploma Thesis; University of West Attica: Athens, Greece, 2023; Available online: https://polynoe.lib.uniwa.gr/xmlui/handle/11400/3968 (accessed on 24 June 2023).

26. Matsoukas, I. Development of an Air Flow Measurement Device Employing a Piezoelectric Membrane and a 3D Printed Structure Containing a Piezoresistor. Diploma Thesis; University of West Attica: Athens, Greece, February 2023; Available online: https://polynoe.lib.uniwa.gr/xmlui/handle/11400/3646 (accessed on 24 June 2023). (In Greek)

27. Matsoukas, I.; Pagonis, D.N. (University of West Attica, Athens, Greece). Personal communication; 2023.