Piezoelectric Power Generation from the Vortex-Induced Vibrations of a Semi-Cylinder Exposed to Water Flow

The aim of this work is to design a piezoelectric power generation system that extracts power from the vibration of a cantilever beam. A semi-cylinder placed in a water stream and attached to the beam is excited into vortex-induced vibrations (VIV), which triggers the piezoelectric deformation. The mechanical system is modelled using parametric equations based on Hamilton’s extended principle for the cantilever beam and the modified Van der Pol model for the bluff body (the semi-cylinder). These equations are simulated using the MATLAB software. The dimensions of the model, the flow velocity and the resistance are treated as design parameters and an optimization study is conducted using MATLAB to determine the combination of optimal values at which maximum power is extracted. The key findings of this research lie in the identification of the effect of changing the design parameters on output power. In addition to the numerical simulation, a finite element analysis is carried out on the bluff body and the hydrodynamic forces and velocity profiles are observed. It is determined that the vibration amplitudes increase with increasing diameter of the bluff body, length of the bluff body and water velocity.


Introduction
Traditional fossil-fuel-based power plants tend to pose many risks to the environment including releasing particulate matter that binds to nitrogen oxides and sulfur dioxide, which in turn causes health problems such as asthma, bronchitis and premature death [1]. In addition, many power plants are built near water streams to be used as coolants. The same water used for cooling is discharged back into the water, only now in its pollutant form, which is extremely harmful to the aquatic life [1]. These risks are only two of many. UCS USA [2] reported that there is a link between running power plants and global warming due to carbon and methane emissions. With the rising pollution rates over the last few decades, there is an urgent need for power generation from renewable sources. In 2017, hydropower was reported to account for 1.8% of Britain's total electricity supply, which made up 18% of Britain's renewable energy supply [3]. Hydropower mainly comes in the form of dams that control water flow and turbines which, in turn, convert stored energy in the water to electrical energy. However, this method of power generation is accompanied by a few risks. Dams have the capacity to impact rivers and the aquatic life by preventing the upstream and downstream migration of aquatic organisms and nutrients and by altering the hydrological flow pattern of rivers [4]. Therefore, a different method of power generation from water flow is examined in this paper. This technique involves the flow of water over a cylinder immersed in a river stream. The flow excites the cylinder into crossflow vibrations, which causes the deformation of a piezoelectric beam attached to it.
Vortex shedding occurs when a flow separates from an elastic bluff body submerged in the fluid [5]. When a cylinder is submerged in a stream of water with velocities exceeding laminar flow, the inertia of the water renders it incapable of negotiating the half of the cylinder facing the downstream. As such, the water flow tends to separate from the top surface of the cylinder, and then peels off in a clockwise motion. This clockwise motion of the fluid is a series of shed vortices, known as the Karman Vortex Street. The same phenomenon occurs at the bottom surface of the cylinder, only the vortices are shed in a counter-clockwise motion [6]. The shedding of these vortices exerts a periodic force on the body, causing it to vibrate transversally with respect to the flow in a simple harmonic motion, provided that the cylinder is restricted such that is does not move along the water flow [6]. These oscillations are known as vortex-induced vibrations (VIV). When the frequency of this motion (the natural frequency) nears the vortex shedding frequency (within ±10%), nonlinear resonance is observed and the amplitude of the cylinder vibrations increases to a peak [7]. This is known as self-excited oscillations VIV. On the other hand, forced-oscillations VIV occurs when the velocity and amplitude of the cylinder vibrations are preset to forcibly produce resonance [7].
In general, recent studies on VIV focus on the effect of certain parameters on the response of a vibrating cylinder and the efficiency of energy extraction. They also investigate different ways to harvest energy from the mechanical vibration of the cylinder.
Xu, Qin, & Gao [8] studied the excitation regions and displacement curves and amplitudes in vortex-induced vibrations of slender, flexible cylinders. Their study shows that there were two excitation regions for different ranges of reduced velocity. The cylinder assumed different displacement shapes in each of the two excitation regions.
Leclercq and de Langre [9] studied the effect of flow-induced bending (in-line deflections) due to the drag caused by VIV on the dynamics of slender cantilever cylinders. It was found that the bending observed in the cylinder exposed to fluid flow switched the system from a single mode lock-in to a multi-frequency response. An impact of this transformation was reduced amplitude vibrations.
Vortex-induced vibrations can be used to generate power. There are different systems that have been built/designed in the past. For example, the cylinder can be mounted on a piezoelectric transducer that converts the deformation caused by the movement of the cylinder into an electric charge [10]. Another method is to attach the cylinder to a linear electromagnetic alternator (a magnet and a coil), where the relative motion between the magnet and the coil induces a current through the coil as per Faraday's Law of electromagnetic induction [10]. There is also the VIVACE (Vortex Induced Vibration Aquatic Clean Energy) system which was patented in 2008 [11]. This transmits the mechanical energy produced by the cylinder to a generator via a gear-belt system. The generator then transforms the mechanical energy into electrical energy [12].
Among the different methods for power generation, piezoelectricity has been proven by Sun et al. [13] and other previous experimentations to produce power in low-velocity flows. The electromagnetic method requires relatively demanding maintenance requirements compared to piezoelectricity [14].
Piezoelectric materials are materials than can produce an electric voltage when exposed to mechanical strain. This is known as the piezoelectric effect. It also works in reverse, that is, the materials tend to deform if a voltage is applied to them [15]. This principle is the basis of piezoelectric transducers. According to Woodford [16], the crystals inside a piezoelectric material are electrically neutral but not symmetrically arranged. When a mechanical stress is applied to it, however, the atoms rearrange due to the asymmetry, which upsets the balance of positive and negative charges and causes an electric potential to appear.
The amount of generated electrical energy is mainly dependent on the strain induced in the piezoelectric layer and on the kind of this material. Different piezoelectric materials can be used to generate power. The materials with the highest piezoelectric properties include PMN-PT and PZN-PT. However, they are highly sensitive to temperature change, susceptible to fatigue and difficult to manufacture. Lead Zicronate Titanate (PZT) is the most commonly used material for piezoelectric energy harvesting since it has sufficient piezoelectric properties while more convenient in terms of sensitivity to temperature, fatigue and manufacturing when compared to PMN-PT and PZN-PT [17].
The PZT ceramics used in energy harvesting [18] are replaced by inorganic and composite materials. Among them, microfiber composite (MFC) material is favoured due to its good bending performance which leads to considerable strains of piezoelectric elements [19]. Shan et al. [20] investigated the feasibility of energy harvesting of an MFC piezoelectric energy harvester (PEH) induced by water vortex shedding.
Modeling and simulation of a piezoelectric generator helps to evaluate the performance of new materials with varying piezoelectric properties. Mokhtari et al. [21] introduced and investigated an electrospun PVDF/LiCl nanogenerator to determine the piezoelectric coupling constant using both experimental and analytical approaches.
Energy harvesting is affected by the orientation of the piezoelectric harvester [34,35]. Coupled damping is the smallest in the vertical arrangement among all configurations so this arrangement produces the largest power.
The most promising method for increasing the output power of a piezoelectric harvesters system is tuning the structural frequency to synchronize with the vortex shedding frequency [36][37][38]. This phenomenon is termed as lock-in. The efficiency of energy extraction through buffeting, galloping or VIV is maximized in the vicinity of lock-in [39]. Moon et al. [40] investigated the effect of cantilever beam mass tuning on the output power and an increase of 508.5% was obtained.
In the other studies different arrangements of the VIV energy harvesting setups have been analyzed. The bluff body can be placed at the upstream side of the piezoelectric in the direction of fluid flow. The influence of wake interference is a defining factor when the harvester is placed behind the bluff body. The lock-in region would not occur if the harvester is arranged in the wake of the bluff body with a small centre-to-centre spacing (1.5 times the diameter). In cases where the spacing ratio is larger than 3, the efficiency would decrease significantly due to the negative effect the spacing would have in terms of VIV suppression. The wake structure is sensitive to the spacing between the harvester and bluff body, the bluff body shape and the flow velocity [41]. The spacing varies over time due to the vigorous streamwise response and as a result the wake interaction varies too. Abdelkefi et al. [42] evaluated the output power of a square piezo aeroelastic energy harvester placed downstream of an oscillating circular cylinder. A wake-induced galloping was observed for the square cylinder at a wide operating range of wind speed.
Suitable geometry significantly improves energy conversion efficiency. It is due to the fact that the highest output power will be obtained when the structural frequency matches the resonant vibration frequency. The VIV amplitude and frequency vary with the cross section. There are several studies which have reviewed VIV for a single circular cylinder [43][44][45][46][47]. The amplitude, frequency, and vorticity distributions are the most important parameters that effectively describe and quantify flow-induced motion (FIM). In particular, the effect of cross section shape on FIM response and energy harvesting of a bluff body was analysed by Lian et al. [48]. This analysis showed that a triangular and a circular cylinder have a higher energy transfer ratio than a semi-cylinder.
Different cross sections have been extensively studied in VIVACE converter or application of VIV generators in wind [49] but results are not applicable for generators with a piezoelectric cantilever beam. Also, parameters may have different effects depending on the selected model. For example, for the VIV model, the onset velocity decreases with the increase of bluff body mass, while for the galloping model, the onset velocity increases with the bluff body mass. The galloping model has two optimum load resistance values at relatively small flow velocities while the VIV model has only one optimum load resistance at all flow velocities.
Beams with other common sections such as square, triangular, and D-shaped have also attracted attention. Ewere et al. [50] developed a model for a square beam and predicted the output power yield. Zhang et al. [51] studied the effect of attack angles on the harvested power of a square beam. It was found that an attack angle of 45 • is the optimal one. In a study of the triangular cross section case [52], a vertex angle of 130 • was recognized as the most preferred. Trapezoidal beams present a greater potential for energy harvesting as compared with a square beam [53]. A D-shaped beam shows a maximum output power of 1.14 mW at a wind velocity of 23.38 m/s [54] whereas peak output power for a T-shaped piezoelectric cantilever is 4 mW in wind speed of 4 m/s [55]. Studies show that other shapes such as three-blade with Y-type cross-section [25] or a fork-shaped cross section bluff body [56] enhance the performance significantly and have a lower onset of wind speed.
The onset galloping speed is sensitive to the cross sectional shape [57]. At a constant resistance load the D-shaped beam achieved the highest harvested power [58]. Also in a study by Sun et al. [13] in a circulating water channel, the maximum power density was achieved through a D-shaped bluff body. Nevertheless, Yang et al. [59] reported that the square beam reached the highest peak power in comparison with the rectangular, triangular, and D-shaped beams.
Previous research that studied the fluctuating cross-flow forces and vortex shedding around bluff bodies with the square beam [60] and a triangular prism [61] experimentally found a significant link between flow orientation and aspect ratio. The decrease in aspect ratio is generally accompanied with increasing vortex shedding frequency.
Sun et al. [13] investigated piezoelectric energy harvesting by VIV and galloping of water whose velocity ranged between 0.2 m/s and 0.54 m/s. In their study, three bluff bodies with different cross-sections were used: circle, semi-circle and isosceles triangle with an 80 • top angle, and the results in terms of vortices formation and power generation were compared. They showed that power output increases with flow velocity and bluff body mass and that the onset velocity increases then decreases as load resistance is increased. Simulations of the experimental setup and the modelling equations showed that the semicylinder was exposed to the highest lift force, which translates to the largest vibration amplitudes, and in turn, the highest capacity for power harvest among the three shapes.
Qureshi et al. [14] studied piezoelectric power generation combined with VIV. The article states that piezoelectric energy harvesters produce power when they are exposed to cyclic pressures in the up and down direction due to fluid flow. They require no maintenance, unless the surrounding conditions vary significantly from the optimal conditions under which they were designed. Piezoelectric energy harvesters are composed of piezoceramic (or piezo-plastic) cantilevers that transform the kinetic energy of the water flow into electrical energy. Their work showed that a 15 mm-long cantilever is capable of producing 0.82 mW only.
The aim of this research is to study the generation of power from the deformation of a piezoelectric transducer mounted directly on the bluff body in a vertical arrangement. The vortex shedding behind the bluff body (semi-cylinder) alters the pressure distribution causing periodic lift forces to act on the bluff body which excites it into crossflow vibrations. On the bluff body, a substrate layer is mounted and attached to a piezoelectric sheet. The piezoelectric material becomes compressed or tensioned and produces electrical energy.
For this study, the system is considered to be placed in a river in Huddersfield, GB. The velocity of the water is approximately 20 cm/s and the river is almost one-meter in depth.
Previous works that discuss piezoelectric power generation induced by vortex-induced vibrations do not optimize the system. Thus, the main objective of this paper is to find the set of design parameters (to be defined later) that allows the extraction of maximum power that the system is capable of producing.

The System
The system consists of a semi-cylindrical bluff body of Silicon material attached to a cantilever beam. Due to its density compared to the density of water, Silicon tends to resist the flow and thus ensures vortex shedding. The cantilever beam is composed of a substrate layer and a piezoelectric layer. The substrate layer is a carbon sheet while the piezoelectric layer is a macro fiber composite (MFC) sheet, which consists of piezoelectric fibers. The top surface of the bluff body is placed below the surface of the water. The cantilever beam is extended into the air. The top surface of the substrate is mounted such that the bluff body and cantilever beam are prevented from moving along the water flow. The flow hits the flat surface of the semi-cylinder and is perpendicular to it. Vibrations are to only take place transversely, with respect to the velocity vector. The cantilever beam and bluff body setup and a dimensioned drawing of the model are shown in Figure 1. The model is also connected to a resistance box. Weirs and notches can be placed behind the setup to increase the velocity of the flow moving towards the bluff body. The most convenient method can be decided upon once the velocity required for power generation is obtained. L tip and D are the length and the diameter of the bluff body respectively. L, L 1 and L 2 are the length of the cantilever beam, the length between the fixed end of the cantilever beam and the starting point of the piezoelectric layer and the length between the fixed end of the cantilever beam and the ending point of the piezoelectric layer. b p and b s are the widths of the piezoelectric layer and substrate layer respectively while h p and h s are the thicknesses of the piezoelectric layer and substrate layer respectively. U if the velocity of the water and R is the value of the resistance specified by the resistance box.

Modelling
The study begins by defining the equations that govern the system. The modelling equations link the vibrations performed by the bluff body and the fluid flow. The equations also link these vibrations to the deformation observed in the piezoelectric layer of the cantilever beam. Lastly, the power extracted from the deformation of the piezoelectric layer is defined. After defining the modelling equations, an optimization study is conducted using MATLAB. Eleven design parameters are identified for this project and the optimization is carried out to determine the set of design parameters that harvests the maximum power. A set of constraints and bounds is also specified for the design parameters.
To find the optimal parameters that would harvest maximum voltage, and in turn power, the model described in The Setup Section of this report is governed by a set of parametric equations based on design parameters. The bluff body has a semi-circular cross-section. These equations are applicable to one semi-cylinder only. The number of semi-cylinders in parallel required to harvest a certain amount of power can then be obtained based on the results of one semi-cylinder.

Mathematical Approach
Using Ohm's Law and the relation between power and voltage, the harvested power can be expressed as where P is the harvested power and V 0 is the harvested voltage. R is the value of the resistance provided by the resistance box. The resistance box is a source of variable resistance. The harvested voltage is defined as [13] where θ p is the piezoelectric coupling term, Ω is the modified frequency of the energy harvester in the modified Van der Pol model, Q 0 is the amplitude of the mode coordinate Q(t) and C p is the capacitance of the MFC layer. Q(t) describes the vibration motion of the bluff body. The capacitance C p is defined as [13] with b p , h p , L 1 and L 2 being the width of the piezoelectric layer, the thickness of the piezoelectric layer, the length between the fixed end of the cantilever beam and the starting point of the piezoelectric layer and the length between the fixed end of the cantilever beam and the ending point of the piezoelectric layer respectively (refer to Figure 1). ε s 33 is the permittivity at constant strain, also known as the dielectric constant.
The electromechanical coupling term θ p , which describes the conversion of mechanical energy into electrical energy, is fully defined by using Equations (4)- (7) shown below [13] where φ is the derivate of the mode shape φ with respect to x, ϑ p is the piezoelectric coupling coefficient, e 31 is the piezoelectric stress constant, d 31 is the strain coefficient of the piezoelectric layer and y is the position of the neutral axis. h s and b s are the thickness and the width of the substrate layer respectively (refer to Figure 1) and E s and E p are the Young's moduli of the substrate and piezoelectric materials respectively. The amplitude of galloping of the bluff body Q 0 is calculated using Equations (8)-(11) [13]. Galloping is the low-amplitude vibrational motion that occurs at low velocities and is incapable of producing vortices that exert negative hydrodynamic forces.
X is the negative linear damping coefficient, w 1 is the first natural frequency, ξ is the mechanical damping coefficient, C is the electrical damping and γ is the cubic hydrodynamic coefficient. k 1 and k 3 are hydrodynamic coefficients and U is the incoming water velocity. The modified frequency is given by [13] with k 2 being another hydrodynamic coefficient. Thus, it has to satisfy the following equation: The hydrodynamic coefficients are defined as shown below [13] where ρ w is the density of water, D and L tip are the diameter and the length of the bluff body respectively and L is the length of the cantilever beam (refer to Figure 1). a 1 is the empirical experimental linear galloping coefficient, while a 3 is the empirical experimental cubic galloping coefficient. C M is the added mass coefficient and in this case, is equal to the density of the bluff body divided by the density of water.
Two main factors affect the frequency of vortex shedding. The diameter of the bluff body is inversely proportional to the vortex shedding frequency, and the flow velocity is directly proportional to the vortex shedding frequency [6]. These relations are depicted by the Strouhal relationship (St = f t D U ) where St is the Strouhal number, f t is the vortex shedding frequency, D is the diameter of the bluff body and U is the velocity of the flow [5]. The relationship between the Strouhal Number (St) and Reynolds Number (Re D ) can be estimated using the following empirical formula [62]: Over a large range of Reynolds Number (250 < Re D < 2 × 10 5 ), the Strouhal Number is almost constant at a value approximately equal to 0.2 [62]. With the Strouhal number being a constant at around 0.2, the vortex shedding frequency can be calculated.
Finding the First Mode Shape at x = L and the First Modal Frequency The modal shapes describe the deformation that occurs in the cantilever beam.
x is the coordinate along the cantilever beam and x = 0 is the point at the top of the beam where it is in contact with the support. The first modal shape can be found using equations Equations (18) where ∅ 11 , ∅ 12 and ∅ 13 are applicable for the conditions 0 ≤ x ≤ L 1 , L 1 ≤ x ≤ L 2 and L 2 ≤ x ≤ L respectively. The coefficients A 11 , B 11 , C 11 , D 11 , A 12 , B 12 , C 12 , D 12 , A 13 , B 13 , C 13 , D 13 , β 11 and β 12 can be found using the following boundary and orthogonality conditions [13]: In the above equations, m 1 = b s ρ s h s and m 2 = b s ρ s h s + b p ρ p h p are the mass per unit length of the substrate layer and the mass per unit length of the substrate and piezoelectric layers respectively, with ρ s and ρ p being the densities of the substrate and piezoelectric layers respectively. The stiffness of the substrate layer EI 1 and the stiffness of the substrate and piezoelectric layers EI 2 are given by the following equations: where y 0 = −y, y 1 = h s − y and y 2 = h s + h p − y. M t is the mass of the bluff body and is the product of the density of the Silicon material that comprises the bluff body and its volume (M t = ρ bb ). L c is the length between the center of gravity of the bluff body and the starting point of the piezoelectric layer (L c = L tip 2 + L − L 2 ). I t is the rotational inertia of the bluff body with respect to the end of the beam. Since inertia is the product of the mass and the square of the length, I t = M t × (L c + L 2 ) 2 where L c +L 2 is the distance between the center of gravity of the semi-cylinder and the end of the beam.

11
EI 1 m 1 The boundary and orthogonality conditions give rise to fourteen complex equations to be solved with the aid of MATLAB. All these equations, along with the cantilever beam mechanical property β 11 = 4 EI 2 EI 1 m 1 m 2 β 12 [13], comprise a system of fifteen equations and fourteen unknowns, to be solved using MATLAB, such that the mode coefficients, and in turn the mode shapes, are obtained. Dimensions, resistance and flow velocity, are treated as design variables.
Equations (1)-(37) used to describe the cantilever beam behavior are based on Hamilton's extended principle, which focuses on the displacement of the system and its time derivatives, as well as the kinetic and potential energies that the mechanical system holds. The equations used to describe the vibrations of the bluff body are based on the modified Van der Pol model which is an oscillatory model with non-linear damping.

Finding the Onset Velocity
The system is to be placed in a river in Huddersfield. The water velocity is low. However, it may be required to increase the water velocity by adding weirs or notches. The onset velocity is the velocity below which no considerable power is produced. The onset velocity is defined as [13] U 0 g = 2(2ξw 1 + C) ρ w Dk g 1 a 1 (38) where k g = ∅ 2 (L)L tip + ∅(L)∅ (L)L 2 tip + 1 3 ∅ 2 (L)L 3 tip . In this paper, the water velocity is treated as a design parameter. However, during optimization of this parameter, the onset velocity if calculated at each iteration to ensure the water velocity used as an input is greater than the onset velocity.

Defining the Constants
The densities of water, the substrate layer, the piezoelectric layer and the bluff body material (Silicon) have the values ρ w = 1 g/cm 3 , ρ s = 8.9 g/cm 3 , ρ p = 5.44 g/cm 3 and ρ bb = 2.329 g/cm 3 respectively. The Young's Moduli of the substrate and the piezoelectric layer have the respective values E s = 110 GPa and E p = 33.336 GPa. Previous experimentation [63] show that the mechanical damping coefficient for a semi-cylinder is ξ = 0.013. The permittivity of free space is ε 0 = 8.85 × 10 −12 F/m. The piezoelectric strain coefficient of MFC (the piezoelectric layer) is d 31 = −170 pC/N. The permittivity at constant strain for the MFC is e 3 ss = 14.04 nF/m [13]. Finally, the hydrodynamic galloping coefficients a 1 and a 3 for the semi-cylinder are equal to 0.05174 and −0.02429 respectively [13].

Numerical Simulation Constraints
There are certain constraints to this design analysis. From the setup's design (Refer to Figure 1), it can be seen that the separate lengths L 1 and L 2 should be less than the length of the cantilever beam. In addition, the width of the piezoelectric layer should be less than or equal to that of the substrate layer. Therefore, the design problem is subject to the following constraints: In addition to the aforementioned constraints, some of the design parameters are also bounded. The substrate and piezoelectric layers can be found in different sizes and thus, the parameters h s , b s , h p and b p are unbounded, except by the aforementioned constraints. The parameter L is unbounded, and the upper limits of the lengths L 1 and L 2 are defined by the value of L. L 1 should be greater than 0 to allow the cantilever beam to be attached to the top support through the substrate layer rather than the piezoelectric layer. The depth of the river environment where this setup is intended to take place is no greater than one meter, which means that the length of the bluff body L tip should be less than 80 cm. The water flows at approximately 20 cm/s. However, weirs and notches can be used to increase the water velocity to around 1 m/s. The diameter D and the resistance R are unbounded. Therefore, the bounds of the design parameters are defined by the following:

Optimization
The design problem is simulated using a number of MATLAB codes. First, the boundary and orthogonality conditions (Equations (21)-(34)) are expanded using a separate code. These expanded equations are then copied into a function code, the aim of which is to find the values of the coefficients A 11 , B 11 , C 11 , D 11 , A 12 , B 12 , C 12 , D 12 , A 13 , B 13 , C 13 , D 13 , β 11 and β 12 . As for the main code, it starts by defining the constants, the parameters that can be calculated from the constants and the design parameters. Then a "for loop" is introduced that simulates Equations (1)- (16), (20) and (35)- (37). Inside the "for loop" is also a statement that calculates the onset velocity (Equation (38)) and checks whether or not this value is greater than the flow velocity. It is observed that when the flow velocity is less than the onset velocity, the harvested power obtained is negligible. Therefore, for simplicity and for a better reading, the code is made to return a 0 value for the power when the onset velocity exceeds the flow velocity.
There are 11 design parameters for this optimization problem: R, D, U, L tip , L, L 1 , L 2 , h s , b s , h p , and b p . Each of these parameters is defined as a vector of values (which adhere to the bounds and constraints) while the other parameters remain constant. The power at each cell of this vector is calculated using a "for loop", and the optimal parameter at which maximum power is observed, is identified accordingly. Several rounds of this operation are performed for the 11 design parameters, until two consecutive rounds produce the same optimal values.

Verification of Mathematical Model and MATLAB Codes
Prior to carrying out optimization, the mathematical model describing the mechanical setup, as well as the corresponding MATLAB codes are verified by comparing the theoretical data produced by the codes with experimental data of literature [30]. The experimental setup used in [30] is identical to the one shown in Figure 1. The values of the design variables used in the experimental setup are shown in Table 1. The MATLAB code is modified to match the eleven design parameters. The voltage V 0 is calculate using Equation (2) Q 2 0 ) and all the subsequent equations that relate it to water velocity. The onset velocity for this experimental setup as per the MATLAB codes is 0.3178 m/s.

Optimal Design Parameters
The maximum power harvested is around 144 mW, with an onset velocity of 49.666 cm/s. The optimal design parameters are shown in Table 2 and Figures 3-13 show the change observed in power over a range of each of the design parameters. The range of each design parameter over which this study is conducted are mentioned later in this section. Some of the graphs are zoomed-in for better plot-readability.
The flow velocity U is varied between 0.2 m/s and 1 m/s with a step of 0.005 m/s. According to Figure 3, the onset velocity is 49.66 cm/s. For flow velocities below this onset velocity, no power is generated. Beyond this velocity, the harvested power increases with flow speed. It can be observed that the maximum harvested power occurs at the upper limit of the flow velocity; that is, at 1 m/s. The resistance R is varied between 500 kΩ and 2000 kΩ with a step of 10 kΩ. Figure 4 shows that for values of the resistance less than 550 kΩ, the model produces negligible power.   For values of R between 550 kΩ and 680 kΩ, the onset velocity is larger than the flow velocity. The maximum power is observed at a resistance of 840 kΩ, beyond which a drop in power in observed. This is also seen in the design equations. The electrical damping C (refer to Equation (9)) increases with the value of the resistance to a maximum and then starts to decrease. The electrical damping then affects the mode coordinate amplitude Q 0 (refer to Equation (8)), which in turn affects the harvested voltage. The source impedance of the piezoelectric layer is defined as Z = 1/ΩC p . This gives a value of Z equal to 843 kΩ, which is almost identical to the load resistance R. This means that the model is in accordance with the maximum power transfer theorem, which states that the highest output power is obtained when the source impedance is equal to the load impedance.
The extracted power is calculated for values of the bluff body diameter D ranging between 40 mm and 60 mm with a step of 0.1 mm. The results observed in Figure 5 for the parameter D show that up to 48.7 mm a negligible power value is obtained, or a large onset velocity is seen. The optimal value of D is 49 mm after which a drop in power output is observed. This is due to the fact that the mass of the bluff body is primarily affected by D. Hence, beyond a certain M tip , the piezoelectric beam starts displaying deformations that produce power less than the optimal power. This means that the hydrodynamic forces of the fluid at the optimal velocity tend to produce different vibration amplitudes in the bluff body as its mass increases. These vibration amplitudes drift farther from the optimal vibration-deformation combination as the mass of the bluff body increases, as seen in Figure 5. The length of the bluff body L tip is increased from 100 mm to 800 mm by 5 mm steps. As can be observed in Figure 6, the model starts to produce positive power for values of L tip being 160 mm or higher, with the maximum observed at 160 mm. When the length of the bluff body exceeds 310 mm, the onset velocity becomes very large compared to the flow velocity and power is no longer produced. Hence, beyond a certain bluff body mass M tip , the piezoelectric beam is rendered incapable of displaying the deformation that produces the optimal power. This is due to the overall inertia of the bluff body reaching a point at this M tip where the resonant frequency is much lower than the forcing frequency caused by the fluid. This renders the vibrations insufficient to produce cantilever beam deformations that lead to the extraction of considerable power. Therefore, a higher flow velocity, and hence larger hydrodynamic forces, is required for increasing values of M tip . The negative power observed in Figure 6 is due to the fact that at certain values of L tip , Equation (12) yields unreal roots, which leads to a negative value of V 2 0 and in turn a negative value for power. The length of the cantilever beam L is changed between 40 mm and 400 mm with 2 mm steps. The same trend observed in the graph of L tip is seen for the parameter L in Figure 7, with the maximum observed at L = 90 mm. It is noteworthy that the mass of the cantilever beam also affects the deformation of the cantilever. In addition, for values of L exceeding 96 mm, the power produced is either negligible or the flow velocity is less than the onset velocity, which means that a higher water velocity is required for the higher mass model, to exert a hydrodynamic force capable of producing bluff body vibrations that cause sufficient piezoelectric deformation. There is noticeable fluctuation in power values as the parameter L 1 (the length between the fixed end of the cantilever beam and the starting point of the piezoelectric layer) is increased from 1 mm to 6 mm by 0.1 mm. This means that the obtained power is highly sensitive to this parameter. A maximum is seen at 4.3 mm in Figure 8. However, the third highest power is observed at 2.5 mm. The difference between the powers at these two values of L 1 is less than 0.0047%. That being said, and given the need to decrease space and budget, the optimal value of L 1 is assumed to be 2.5 mm. With the values of L and L 1 being constant, changing L 2 (the length between the fixed end of the cantilever beam and the ending point of the piezoelectric layer) means changing the length of the piezoelectric layer. L 2 is varied between 10 mm and 100 mm with a step of 1 mm. Figure 9 shows that the water velocity becomes sufficient at L 2 = 63 mm and maximum power is obtained at L 2 = 70 mm. For values greater than this, the power produced decreases and fluctuations are observed. This drop is due not only to the mass of the model, but also the piezoelectric material's maximum electricity-generating capacity and the effect of L 2 on piezoelectric deformation. The width of the piezoelectric layer b p is increased from 10 mm to 50 mm by 0.5 mm steps and the power extracted is reported at each step. For values of b p less than or equal to 24 mm, the water velocity is insufficient to produce considerable power, as observed in Figure 10. The optimal value for b p is 31.5 mm. For a value of b p = 29 mm, there is a significant drop in harvested power. Investigating the MATLAB code shows that the drop is due to a decrease in the electromechanical coupling term θ p . This term depicts the conversion ratio between mechanical energy and electrical energy. In turn, this drop in θ p is due to a low value of the term φ (L 2 ), which is the derivative of the mode shape of the beam with respect to x at a length L 2 . It can be concluded that at b p = 29 mm, and at the optimal piezoelectric length and thickness, the deformation in the beam does not cause the electric charges to align such that an electrical output is produced, hence the low value observed for θ p . The thickness of the piezoelectric layer is varied between 0.1 mm and 10 mm with a step of 0.1 mm. According to Figure 11, the maximum energy harvested occurs at an h p value of 0.3 mm, and the critical dimension above which power becomes negligible due to onset velocity is 0.4 mm. It is noteworthy that power is extremely sensitive to this parameter as a 0.1 mm change causes the power to drop from maximum to minimum. This analysis leads to the conclusion that a very specific thickness is required for the electric charges in piezoelectric materials to properly rearrange upon the given deformation such that they produce considerable power.  Figure 12 shows that the optimal value of the width of the substrate layer b s is 34 mm. This optimality is obtained for the parameter b s ranging between 31.5 mm (the optimal width of the piezoelectric layer) and 50 mm, studied every 0.5 mm. Values of this dimension below 34 mm, produce positive power which increases up to the maximum at b s = 34 mm. For values of b s of 37 mm and higher, power output is negligible. This is because the width of the substrate affects the bending capacity of the cantilever beam and thus, the deformation of the piezoelectric layer. The optimal value of the thickness of the substrate layer h s , ranging between 0.1 mm and 10 mm with 0.1 mm steps, is observed in Figure 13 to be 6 mm. At 6.4 mm, and for a velocity of 1 m/s, the power output is reduced to zero, due to the increased thickness reducing the bending capacity of the beam.  Figure 14 shows the 3D figure of the optimal model and its dimensioned drawing. The mass of the entire model is 518.23 g and its volume is 169, 857.15 mm 3 or 169.857 cm 3 . These values are obtained from SolidWorks and exclude the top support from which the setup hangs. Using Equation (17), the vortex shedding frequency is found to be 4.08 vortices/second.

Optimal Model Specifications
Using the values of mode coordinates and mode shapes extracted from the code, the maximum tensile strain is found to be 0.015 which is less than the recommended maximum tensile strain (0.95 on average). This means that the solution to the optimization problem is acceptable. The mechanical power MP produced by the vibrations of the bluff body is given by [64]: where C L is the lift coefficient and equal to 0.5 based on experimentation on setups with similar Reynolds number [64] and w s is the vortex shedding frequency. This formula gives a mechanical power of 1.24 W or MP = 1240 mW. With an electrical output of 144 mW, the electromechanical efficiency is found to be 11.6%. Table 3 shows the densities and piezoelectric properties of different piezoelectric properties. The materials include the MFC defined in this study. It also includes PZT-5A, PZT-5H and PMN-PT (30% PT). The properties of the three new materials are extracted from [65]. The results of this comparison in terms of power output and onset velocity are shown in Figure 15.

Comparing the Performance of Different Piezoelectric Materials in the Optimal Setup
When comparing the performance of the different materials, Figure 15 shows that the MFC produces the highest power, followed by PMN-PT (30% PT), PZT-5H and then PZT-5A. However, in low-velocity flow environments, PZT-5H is the most convenient since it has a low onset velocity. Among the four materials, PZT-5A has the least performance in terms of power output and onset velocity.

Computational Fluid Dynamics (CFD)
The amplitude of vibration of the bluff body depends mainly on the hydrodynamic forces exerted by the fluid on the body. To determine hydrodynamic forces, a CFD analysis is conducted using ANSYS Fluent. The inlet velocity is specified, and the solution is calculated at 0.001 s per time step for 8000 time steps.

Mesh Convergence
A mesh convergence study is conducted to ensure that the elements of the mesh are small enough to accurately compute results without being needlessly small. The results show that the lift force has changed less than 2% for the element size of 0.04 m and higher ( Figure 16).

Domain Convergence
The domain convergence study is conducted to determine the smallest domain dimensions that produce accurate results which prevents the turbulence at the walls of the domain interferes with the actual solution. The domain convergence is performed based on the lift force. The convergent width is first identified by increasing the width gradually until the lift force converges. The same process is repeated for the length at the convergent width. The results are shown in Figures 17 and 18.  The difference in the lift force between this dimension and all subsequent dimensions is less than 2%. In addition, this dimension is more than 20 times the diameter of the bluff body, which ensures that it is large enough in width and length for the flow to normalize, for turbulence interference to be avoided and for the software to accurately compute results.

CFD Results
CFD analysis is conducted on the optimal model and three models. Table 4 summarizes the parameters of the different models. The variation of the lift and side forces exerted on the bluff body over time are plotted in Figures 19 and 20 respectively for the four models. Figure 21 shows the velocity contours that allow for observing the vortex shedding occurring in each model.
There is a common trend among the plots for all models. The lift and the side forces fluctuate between positive and negative values for all models, and this is due to the von Karman effect (Figures 19 and 20). The lowest side and lift forces are observed in Model C where the flow velocity is half the optimal value and D and L tip are at their optimal value.
Negative cross-sectional area of the semi-cylinder which is affected by both D and L tip . On the other hand, the lift force is exerted on the horizontal cross-sectional area of the bluff body which is affected by variations in D only. Thus, theoretically, the highest side force would be expected to be seen in Model A or Model B, while the highest lift force would be expected to be observed in Model A. The simulation results, however, are not in complete agreement with these theories.
In general, the highest lift force is observed for Model A, followed by Model B and then the Optimal Model, and the most uniform fluctuation is observed in Model A. Similarly, the highest side force is observed in Model A, followed by Model B and then the Optimal Model. These comparisons, however, are not true for all points in time.
The reason for the discrepancy between the theoretical and the simulation results is the fact that the theoretical formulae for hydrodynamic forces are independent of time. On the other hand, the simulation results shown in Figures 19 and 20 are the products of transient analyses which change with time and are highly affected by the turbulence caused by vortex shedding behind and in the wake of the semi-cylinder. Changes in velocity are also observed around the semi-cylinder during vortex shedding which also causes variations in the hydrodynamic forces.
The velocity contours are shown at a specific flow time of 1 s in Figure 21. Model C shows the smallest vortices, due to its low lift and side forces. Model A, Model B and the Optimal Model show vortices similar in size.
These results may seem a bit contradictory with the numerical MATLAB optimization results. While it may seem that Model A produces the largest vibration amplitudes, this does not mean that it is the optimal model. This is because the optimal model is the one that produces the most power and not the one that has the highest vibration amplitudes. Power generation is dependent on the piezoelectricity of the setup. There exists a certain deformation of the cantilever beam, caused by the bluff body vibrations, which produces maximum output power. This deformation is not taken into consideration in the computational fluid dynamics analysis.

Conclusions
In this research, a power-generating system is designed. The power generation is the outcome of the deformation of a piezoelectric beam. This deformation is caused by the vibrations of a semi-cylindrical bluff body immersed in a river and exposed to water flow.
The design methodology begins by defining the modelling equations, which are functions of eleven design parameters. This mathematical model is outlined in MATLAB codes, which in turn are validated by experimental data.
An optimization study is conducted for the design parameters and the maximum harvested power is found to be around 144 mW. The values of the optimal design parameters where maximum power is extracted are identified.
The key findings of this research lie in the identification of the effect of changing the design parameters on output power. The plots of the design parameters L tip , L, b p and h p versus power tend to have sharp peaks or drops around the optimal value, meaning that accuracy is required for these parameters. Other parameters such as R, D, L 2 , b s and h s have a smoother variation curve and a less drastic effect on power around their optimal values. L 1 fluctuates unstably throughout the variation plots. Lastly, U has a lower threshold for producing power, known as the onset velocity, above which power tends to increase almost linearly, which proves that increasing this parameter is most effective in increasing output power.
In addition to the numerical simulation, a finite element analysis is carried out on the bluff body and the hydrodynamic forces and velocity profiles are observed. It is determined that the lift and side forces, and hence, the vibration amplitudes increase with increasing diameter of the bluff body, length of the bluff body and water velocity. This also translates into larger vortex size.
The future of this project lies in conducting an investigation to check the possibility of replacing the piezoelectric beam with a hybrid piezo-electromagnetic beam to produce more power. Furthermore, a real-life implementation strategy must be put in place, with a feasibility study of implementing the setup in high velocity environments (such as oceans) to produce considerably higher power outputs.