Introduction
As stated by the International Energy Agency 1, energetic consumption has increased considerably in the last decades, and significant increases are estimated for the next years. Nonrenewable resources are the main source of energy generation in countries that do not have the necessary technologies to take advantage of other less polluting energy sources.
This reality encourages the development of clean and accessible technologies to use renewable energy sources 2-5. Harvesting energy from vibrations is an alternative to taking advantage of non-conventional energy sources, which have been gaining space in the last decades.
Vibrations as a source of energy can be harnessed by conversion mechanisms such as electromagnetic, electrostatic, and piezoelectric transduction 6-9. Piezoelectric transduction is consolidating as one of the most promising methods in electric energy micro-generation, due to its high power density; moreover, it does not show any electromagnetic interference 10,11.
The studies on piezoelectric materials have allowed the production of composite materials with excellent properties, focused on obtaining the maximum equipment or systems performance, which has led these materials to have a great variety of applications in different fields of science 12-14.
Wind’s kinetic energy is able to produce induced vibrations of vortex, slipstream, flapping, and gallop types 15. In the scientific advances in the use of aero-elastic vibrations and the development of several prototypes, the L-type 16, circular cylinders 17, and triangular leaf shape 18 stand out. The connection in parallel of several harvesting models 19 has also been experienced, among others.
There exists a remarkable variety of harvesting devices models that use basic geometric forms (square, circular, triangular) as surfaces of interaction or of kinetic energy capture, but the use of geometries of certain complexity in panels focused on rising the energy or vibrations harvesting can have a positive effect on the power output; this fact increases performance and reduces in a certain way, the economic costs of production. The selection of profiles is made using criteria such as design simplicity, economic factors, time for manufacturing, size, and maximum stresses, or pressures, among others.
The high variability of the vibration frequencies can affect the harvesters’ performance; therefore the study of these scenarios is of fundamental importance in order to estimate the average electric energy production 19.
Technologies and devices of micro-consumption that operate in daily activities, point to the development of energy harvesters, considering these as one of the main sources of energy for future smart cities, where the use of endogen resources is a fundamental aspect.
The advances in computational tools for design assisted by computers, allow studying physical processes of high complexity. An example of this fact is computational fluid dynamics (CFD), which is capable of analyzing with high precision, the effects of a fluid on an element considered as a frontier or obstacle in the computational domain, independently of its geometry 20-23.
Along with these computational advances, the development reached by numerical methods has allowed the coupling of fluid problems (CFD) with modal analysis (FEM), making possible the identification of zones of higher stresses due to fluid pressure. It also allows the analysis of frequencies in order to investigate the effects of the different vibration modes of a particular element or system.
Most of the CFD models developed on energy harvesters use different flow velocities, being frequent some values within the range of 1.8 to 23 m/s. The geometrical characteristics of the models are different between authors; still, the obtained power outputs varying from 0.001 mW to 4 mW can be related, which are dependent on the flow velocity, the transduction method, and the geometrical parameters used 24,25. CFD results are often not used as a basis for further studies with other techniques such as modal analysis; consequently, in general, the displacements and oscillation frequency of the harvesters are the result of experimental validations 26,27.
Based on the above-mentioned background, this work aims to the simulation of the mechanical response of a model of energy harvester by wind loads, using the method of computational fluid dynamics (CFD) coupled with the modal analysis (FEM).
Materials and Methods
2.1 Description of the phenomenon to be modeled
To take advantage of aero-elastic energy vibrations from a uniform, horizontal air flux, a mathematical model CFD/FEM is developed for a prototype of an energy harvester from the vibrations generated by wind loads.
The harvester’s organ response in function of the different vibrational modes (modal analysis) in FEM is simulated. Subsequently, an analysis of the virtual model mechanical response is developed to define its resistance to loads or fluid pressure CFD and, at the same time, to identify zones of maximum stress concentrations, strains and displacements where the piezoelectric material could be placed and generate an output voltage and usable power for the energy supply of devices or sensors.
2.2 Virtual Model of the harvester organ
The object of study is a piezoelectric energy harvester, consisting of a panel in the form of a wave (wave-shaped by the wind in a fluid) and a rectangular piezoelectric patch. The patch resides on the panel’s internal surface (Figure 1). The organ’s virtual model was generated using the SolidWorks® 2020 software, by some utilities, including sketch, equidistance, and solid operations.
2.3 Definition of the harvester material properties
Typically, in a real scenario, the energy harvesters are in corrosive and abrasive environments, a fact that conditions the selection of materials for their construction, besides the economic cost that the material can represent when having better mechanical properties to give a longer service life. Table 1 shows the materials’ physical and mechanical properties required by the models for the simulation of the harvester prototype functioning and operation.
The piezoelectric material voltage output is conditioned by its properties and the applied force, in piezoelectric materials of a disk with diameters higher than their thickness, in a proportion of 5:129. The theoretical or calculated output voltage can be obtained by means of Equations 1 and 2:

Where t is the thickness, m; d is the diameter, m; g33 is the piezoelectric voltage constant, Vm /N; F is the force on N;V PD is the output peak voltage in V.
Table 1 Physical and mechanical properties of the semi-curved panel and piezoelectric material
| Parts | Properties | Value |
|---|---|---|
| Semi-curved panel (Copper) | Elastic module (MPa) | 200 |
| Poisson Coefficient | 0.35 | |
| Density (kg/m³) | 8,9 | |
| Tensile or traction limit (MPa) | 413 | |
| Elastic limit (MPa) | 172 | |
| Thickness (mm) | 0.2 | |
| Piezoelectric material (PZT NCE51) | EDensity (kg/m³) | 0.0078 |
| Piezoelectric voltage (Vm/N) | 0.0027 | |
| Diameter (mm) | 30 | |
| Thickness (mm) | 1 |
Equation (2) expresses the voltage as a function of the applied force. With the calculation of V PD we proceed to calculate the efficient voltage, as shown in Equation 3:
Where V Ef is the efficient or effective voltage in V. Finally, it is possible to know the electric power from the fundamental equations in Ohm’s Law.
2.4 Model in CFD for the simulation of the effects of wind loads
The simulation made by CFD will allow investigating the factors influencing the energy harvest produced by the wind loads, such as the wind velocity variation, the aerodynamic pressure, shear force, turbulence, and environmental conditions. The CFD model was developed in the SolidWorks® software using the complement Flow Simulation.
A study of external type fluid is defined, where the computational domain is adjusted to the panel’s geometry; the interactions between the structure and the fluid (velocity variations), and the surface objectives for the study of the stresses are established (Figure 2).
The computational domain considers that the laminar type fluid with kinetic energy moves at a constant velocity in the direction -z, with an ambient temperature of 298 K and a pressure of 101.325 kPa, which is within the average ranges for the city of Portoviejo, Ecuador, a zone where the study is performed 29-31.
Concerning the wind velocity, a variation range of 3 to 21 m/s was established to represent this variable’s variations in the different zones of the Republic of Ecuador 30,32,33.
Model’s discretization
The model’s discretization is made by the Cartesian mesh method of the Solid Works Flow Simulation, which gives great advantages in generating the mesh algorithm considering the original CAD model data (virtual model).
The method generates three types of cells that represent the solids, the fluid, and partial cells that represent the intersection with the immersed solid. These cells are rectangular divisions (cuboids) of the computational domain. For the meshing, a manual configuration was used for both the virtual model and the fluid, with an advanced channel refinement; the minimum size of space was 1 mm and a radius factor of 4, resulting in a total of 243 607 cells for the fluid, 15 542 cells of the solid and 20 641 cells in contact with the solid. Figure 3 shows the mesh calibration process, and a number of elements are chosen within the range of 17 200 to 21 200, because the results are independent of the mesh.
2.5 Coupling of modal analysis by the finite elements method (FEM)
In order to determine the harvesting response in the different vibration modes of the harvester work organ, a modal analysis was developed in which there were also determined the model’s natural frequencies and the contributions of the organ mass, in each one of the coordinate axes, where the deformation mode can define the most adequate scenario for the energy capture from the wind loads. From the interaction structure-fluid results, the four lower modal forms were analyzed.
The computational tool allows the coupling of the results of a model CFD with a modal analysis in FEM to link aerodynamic pressure loads. The general procedure for performing a CFD/FEM study in SolidWorks®, is shown in Figure 4.
2.6 Harvester’s physical integrity determination
With the fluid pressure and shear stress calculations in three dimensions (CFD), we proceed to study the mechanical response of the harvester panel’s material, produced by the
action of the mechanical stress simulated by the SolidWorks Simulation software component. By means of this method application, the harvester’s physical integrity against the different pressure magnitudes is evaluated. For the frequency studies, the lower ends of the panel are fixed (Figure 5), allowing the panel to oscillate in the function of the applied loads. Finally, the model discretization was developed, starting from a mesh based on a curvature composed of quadratic elements of a high order, which generates a total of 84 709 nodes and 50 175 elements.
Results and Discussion
3.1 Results of the organ-fluid interaction (CFD)
The harvester organ’s geometry allows the opposition to the kinetic energy transmitted by the wind load, generating pressure on the panel’s surface in the central zone, that reaches a maximum of 101 716 Pa and of 101 336 Pa, for the wind velocities of 21 and 4 m/s respectively, in the zone of direct impact with the wind (Figure 6a).
The wind velocity experiences a small drop of 2 m/s on the surface of the direct impact of the wind on the panel, but in the superior part of the structure (vortex), a small velocity increment is observed (23.48 m/s), while that for the velocity of 4 m/s, it is of 5.859 m/s. Both increments are a result of the fluid accumulation and the low resistance that the organ’s geometry generates in the zone. (Figure 6b).
The relative pressure shows negative values of -244.62 Pa in the panel’s superior zone; that is, a small vacuum is generated, which initiates vorticity; also, in the impact zone, we observe how the pressure rises to 391.75 Pa. This trend is similar to the velocity of 4 m/s, having values of -10.87 Pa in the vortex zone and 11.51 Pa in the impact zone (Figure 6c). Negative values indicate that there is a deflection causing traction, resulting in an expansion of the material, while positive values represent compression or contraction caused by the wind load. On the other hand, the shear stress of the fluid action on the surface is in the range of 2.44 to 4.87 Pa for the maximum velocity, while for the minimum velocity, the stress is 1.55 Pa (Fig 6d). The friction coefficient decreases as the velocity increases, in this way, for the maximum velocity, the friction coefficient is 2.26 for the maximum velocity, and at 4 m/s, the coefficient is 6.47.
Additionally, pressure, shear stress, and relative pressure variables are dependent on the velocity variation; the higher the magnitude of velocity, the higher the resultant of the other variables.
Since the harvesters left and right surfaces have different geometries, the pressures and the other variables can vary depending on the section with which the fluid interacts. The analysis of the relations between dynamic pressure and wind velocity, considering the direction from which the wind blows, allowed obtaining the models that relate both variables (Figure 7). These models evidence the existence of a third- order polynomic relation in all the analyzed cases, showing that, as the fluid velocity raises, the dynamic pressure exerted by the fluid increases in a non-linear way. In cases 1 and 3, it is observed that for velocities higher than 12 m/s, the non-linear relation between both variables is accentuated. The correlation coefficient reached high values of R2 = 0.9926, 0.989, and 0.9874, and the errors of 0.051, 0.029 and 0.034, respectively for each one of the cases.
The maximum pressures were obtained when the wind’s direction acts in the negative direction of the z-axis. The simulation shows that the panel geometry can accumulate a higher aerodynamic pressure with respect to the curved panel 34.
In Figures 7b and 7c, it can be clearly appreciated how the capture geometry (panel) opposes the laminar flux and modifies its displacement that, in principle, is made in an ordered way with a homogeneous velocity between fluid layers. The interaction produces a decrease in the pressure on the panel backside, modifying the fluid’s direction and velocity, forming a forced vortex in a clockwise direction. The variations and location of pressures show significant relations with other models of harvesters that also take advantage of the wind loads to oscillate 17,22,35. This vortex concentration (turbulent flows) after the geometry, and allows us to say that, part of the air’s kinetic energy is transformed into aerodynamic pressure; that is, the flux changes to turbulent below the panel maximum height.
The resulting pressure of the interaction between the fluid and the harvester organ produces small displacements in the geometry, that compress and extend certain zones of the panel, producing vibrations that can be used later.

Figure 6 CFD simulation results, velocity of 21 m/s. (a) aerodynamic pressure, (b) velocity, (c) relative pressure, and (d) shear stress
3.2 Modal Analysis FEM results
Results show the natural frequencies, periods, and amplitudes of the vibrations when the harvester organ is impacted by the wind load (Table 2). The three first resonance modes have low natural frequencies compared to the fourth mode, whose oscillation period is shorter. With respect to the oscillations’ amplitudes, it can be observed that the second modal form, has major displacements (Figure 8).
To get the most out of displacements, it is necessary to identify the higher amplitudes in the axis “y” and “z”, because this, starting from the wind pressures, generates flection and torsion stress that the piezoelectric material can transform into electric output. In accordance with this fact, the more feasible vibration modes in low frequency are the first, the second, and the third modal forms, while the fourth modal form can also be used, but has a higher frequency.
Some design considerations for these harvesters recommend avoiding excessive torsional stresses, which could produce significant deformations, and the frequency of excitation in the harvester should be close to or coincide with the natural frequency in order to reach the resonance state that generates an increase in the amplitude of the oscillation and the electric generation.
The natural frequency values are higher than the ones found by 34; however, because the model is in a dynamic environment, it was considered to link the pressures generated by the fluid, which, in a certain way, induce displacements, deformations, and tensions that modify the modality characteristics, increasing them.
The appropriate zones to place the piezoelectric transductor material are those where the vibration amplitudes are higher; that is, in the zones where the maximum amplitudes are reached. It must be remarked that for each one of the natural frequencies, the location of the zones where the maximum amplitudes are reached, changes.

Figure 8 Vibration modes with an overlap of the model on the deformation, (a) first, (b) second, (c) third, and (d) fourth
Table 2 Frequency modes parameters
| Frequency mode | Natural frequency (Hz) | Period (s) | Amplitude (mm) | |||
|---|---|---|---|---|---|---|
| (x-axis) | (y-axis) | (z-axis) | Resultant | |||
| 1 | 69.013 | 0.01449 | 0.0037 | 0.915 | 1.178 | 15.51 |
| 2 | 88.872 | 0.011252 | 13.28 | 13.01 | 12.31 | 19.41 |
| 3 | 99.482 | 0.010052 | 0.022 | 1.487 | 13.40 | 14.42 |
| 4 | 253.54 | 0.0039441 | 0.021 | 11.52 | 12.22 | 14.11 |

Figure 9 Mass participation factor associated with the excitation frequencies and possible piezoelectric material location for first (a), second(b), and third (c) modal forms with respect to the z-axis
Figure 9 shows the percentage of cumulative mass participation (CEMPF) for the different frequency modes. For the x-axis, the three last modal forms have the same CEMPF of 23.5%. With respect to the y-axis, the CEMPF values for the four modal forms are 42%, 65%, 72%, and 77%, respectively. In the z-axis, the higher participation percentages are obtained for the frequency modes; the third and fourth modal forms stand out, where the CEMPF is 124% and 130%, respectively. This result shows that the piezoelectric material location must be done along the z-axis or the harvester organ’s longitudinal axis, in the zones where the major displacements are observed.
3.3 Results of the structural integrity check of the element or organ
The stress distributions (Figure 10a) show that the maximum von Mises stresses are located in the central zone of the panel, where the wind load impacts directly, as well as in the model’s posterior support. These displacements are in the range of 0.015 to 0.018 mm (Figure 10b). The maximum displacements are able to produce significant flection stress on the piezoelectric material, if its location is near the model’s center.
The unitary deformations in the structure allowed showing that the deformations caused by the wind load are in the range of 0.000000941 a 0.00000219 (Figure 11a), so the panel’s structural integrity is not compromised. These deformations are located directly in the zone of wind impacts with the panel’s surface, which corroborates or indicates the zone where the piezoelectric material should be placed.
A security factor distribution analysis shows that it reached values Fs=44, which evidences that the harvester resists the loads as well as the deformations that arise from wind loads, without affecting its structural integrity. Despite the harvester having these resistance reserves, we decided to maintain the panel’s dimensions because of the significant contributions to the three-axis displacement that a higher mass generates, besides giving a significant wind impact area.

Figure 10 (a) von Mises stresses distribution and (b) displacement, air flux z-axis negative direction
3.4 Mechanical stress transduction
In figure 12, it can be observed that, as the wind load applied on the panel increases, the peak voltage and the efficient voltage of the PZT NCE51 increase linearly, showing a directly proportional relation. For the case of power, a non-linear relation between applied force and efficient voltage is shown, evidencing that as the load rises, the power increases in a polynomic manner, as a result of the vibrations generated in the harvester’s panel. The correlation coefficient reached values of R2=1 for the analyzed cases. It is evident that the location of the piezoelectric material in a different section, of higher stresses, produces higher values of voltage output and power. The use of piezoelectric PZT exhibits higher output power and efficiency, as mentioned by several authors 34-36.
4. Conclusions
The coupling of the results of CFD simulation with the computational simulation by FEM, allowed simulating the wind pressure effects on the harvester’s work organ, as also to know the vibration modes most favorable for the electricity generation, in a range of wind velocities of 3 to 21 m/s.
The CFD simulation results allowed identifying the zones where the pressures concentrate, which vary from -0.244 to 0.391 kPa, as a product of the load impact imposed by the action of the wind on the panel’s surface, as well as to determine the maximum deformations.
Maximum deformations in the three analyzed axes are originated in the second mode of vibration, which makes it the most adequate for electricity generation in the simulated operating conditions.
The piezoelectric material must be placed in the direction of wind action, because this is the direction where the panel suffers the maximum deformations.
The harvester is able to preserve its structural integrity during its operation in the simulated wind load regime.
























