Broadband vibration energy harvesting from a non-deterministic system: Performance of different piezoelectric patch shapes

Harvesting energy from ambient structural vibration using piezoelectric materials gained massive interest in the past decade. Piezoelectric harvesters can be incorporated in many applications; however, one of the main challenges to become widely adopted is to optimize their design for maximum energy harvesting. In this paper, we investigated energy harvesting from a piezoelectric patch that is attached to a non-deterministic thin plate vibrating in bending. Energy harvesting from six patch shapes (differing in the number of edges) was examined through a coupled-field finite element model. The thin plate was simply supported with nominal geometry and material properties. The plate’s dynamics were randomized by randomly distributing point masses on its bottom surface; this made the plate a non-deterministic subsystem. The design optimization was performed by changing the shape of the piezoelectric patch and analyzing the ensemble response of the electrical potential across the piezoelectric patch. The results show that piezoelectric patches with an even number of edges exhibit higher performance in terms of energy harvesting.


Introduction
Modeling the vibration of structures in the high-frequency range is challenging. On the one hand, analytical approaches have their limitations, and on the other hand, the computational cost of finite element models can become prohibitively expensive (especially at higher frequencies). Thus, using deterministic methods becomes less favorable and energy-based methods, such as Statistical Energy Analysis (SEA), come to the forefront [1]. In particular, energy-based methods such as the SEA can be used when the modal overlap factor (MOF) exceeds unity. The MOF is defined as the ratio between the modal bandwidth at the half-amplitude and the average modal spacing [2] and it can be calculated as: where n is the modal density, i.e., the number of modes per unit frequency, h is the modal loss factor, and w is the natural frequency in rad/s. For plates, the modal density n can be expressed analytically as: where S is the surface area, and r, h, and D are the density, thickness, and flexural rigidity, respectively. The flexural rigidity is expressed as: where E and u are the modulus of elasticity and Poisson ratio (respectively) of the plate's material. The MOF can be used to establish three frequency ranges: low-frequency range ( < MOF 1), mid-frequency range ( < < MOF 1 2), and high-frequency range. In the low-frequency range, individual modes are wellspaced, and the response exhibits distinct peaks. The dynamic response of the structure in the low-frequency Any further distribution of this work must maintain attribution to the author(s) and the title of the work, journal citation and DOI. range is not sensitive to uncertainties; hence, such a response can be modeled with element-based methods such as the finite element method (FEM). In the high-frequency range, the contribution of individual models combines in the response and resonant peaks are no longer distinct. Figure 1 shows a representative frequency response function (FRF) that exhibits the three frequency ranges. In the mid-frequency range, the modal overlap and/or loss factor are higher. This intermediate frequency range holds the characteristics of low-frequency; however, the modal response begins to overlap, and it is not clear which modes combine. At the mid-frequency range, a mix of dynamic behavior is encountered where some sub-systems are large compared with a wavelength and thus well modeled by SEA, while other sub-systems are short compared with a wavelength and are not well modeled by SEA [3].
In the low-frequency range, low modal overlapping occurs, the modes are distinct and well-spaced; hence, the system response has high variance. As the MOF increases, the variance of the response becomes very small as high modal overlapping occurs and no single mode dominates the system response. Consequently, structures vibrating in the high-frequency range are typically referred to as non-deterministic substructures (Non-DS). These structures produce a vibration response that is sensitive to uncertainties and can be analyzed using SEA [4]. The real-life response found in Non-DS can be achieved by randomly distributing point masses and/or point springs on the structure to introduce uncertainties to the response. Uncertainties can also be simulated by means of variabilities in structural damping, thickness, material density, and elasticity [5]. The utilization of FE and SEA for hybrid vibration response analysis has been incorporated in many researches in an attempt to model high-frequency vibrations. For example, SEA and FE analyses (FEA) were used to obtain the ensemble-average of the time-averaged vibrational response of dynamic systems [6]. The FEA-SEA method was later extended to obtain the structural responses of systems consisting of FE deterministic components and SEA components under impulsive and transient loadings [7]. FEA and SEA were also proposed as a hybrid technique to optimize the accuracy of response prediction of beam-plate structures in the mid-frequency range [8].
In the past decade, research into using piezoelectric energy harvesters (PEHs) to convert mechanical energy into electrical energy has attracted considerable interest from researchers. PEHs exhibit high energy conversion, are low cost, and can be easily fabricated. PEHs have been incorporated in various applications, including civil structures such as buildings and bridges [9], aerospace [10], medical devices [11], and even in human body movements [12]. Numerous investigations have been conducted with the aim of optimizing the performance of PEHs; this includes optimizing the mechanical and electric components. Shafer et al [13] considered the piezoelectric bimorph energy harvesters in beam applications. They optimized the thickness of the piezoelectric layers relative to the total beam thickness, thus achieving more power output. Chavez et al [14] optimized the performance of PEHs by thermal-mechanical coupling. Deng et al [15] investigated the performance of a cantilevered type PEHs and optimized the damping ratio and electromechanical coefficient for maximum power output. Also, Song et al [16] studied a cantilevered PEH and proposed an optimization strategy to achieve the maximum power output by adjusting the piezoelectric properties and the geometry of the piezoelectric layers. Fan et al [17] proposed a hybrid energy harvester that maximizes the power output by integrating two piezoelectric cantilever beams and an electromagnetic energy harvester to scavenge energy from low-frequency and bi-directional excitations. Another optimization strategy was proposed by Kim et al [18] to optimize the performance of a cantilever-based PEH by devising a method to prevent voltage cancellation and ensure that the system operates in a wide frequency band. The cantilever-based energy harvester in Kim et al [18] operates by shifting its resonant frequency. This was achieved by adjusting the location of the PEH based on the expected variations to the phase angle of the cantilevers. These researches have explored the performance optimization of cantilevered type PEH at the dominant resonant modes that are evident at the low-frequency ranges. However, the performance of PEH in the high-frequency range where the response is more challenging to model due to response uncertainties has not been investigated.
In this article, the performance of the ceramic type piezoelectric energy harvester attached to a Non-DS thin rectangular plate was investigated using FEA. The performance of the PEH attached to the Non-DS plate was investigated by employing an ensemble average response. Piezoelectric patches (which can have any of six shapes differing in the number of edges) were attached to the simply supported thin plate. Randomly located point masses are distributed on the plate to introduce uncertainty, making the system a Non-DS. The system was subjected to a point force excitation, and the electric potential across the patch was obtained. Later on, the ensemble responses were analyzed to study the effect of piezoelectric shape on the developed voltage.
This paper is organized as follows. Section 2 presents the construction methodology of the finite element model of the piezoelectric energy harvesters attached on a non-deterministic rectangular plate including, modeling, material properties, meshing, and boundary conditions. As well as, a description of the piezoelectric effect equations adopted in the study. Section 3 compiles the results obtained from the simulation studies and presents the analysis performed on MATLAB to investigate the response behavior. The paper ends with a conclusion that includes a summary of the work performed throughout the study and presents the findings of the research.

Methodology
ANSYS Mechanical APDL ® (ANSYS ® ) was used to simulate the dynamics of the thin plate and the piezoelectric patches attached to it.

Piezoelectric model
The PEH utilizes the piezoelectric effect to generate an electric potential in response to applied mechanical strain on the plate's upper surface, as illustrated in figure 2. The direct piezoelectric effect is used to describe the interaction between stress, strain, charge, and electric displacement. The piezoelectric effect can be described in strain-charge form using: In equations (4) and (5) The direction of positive polarization coincides with the z-axis of a rectangular system of x, y, and z axes (as shown in figure 3). In the compression mode, the coefficient d 33 would indicate the polarization generated in the z-direction per unit of mechanical compression stress { } T applied in the z-direction to the piezoelectric body, or the induced strain in z-direction per unit electric field applied in the z-direction. The transverse mode involves the coefficient d 31 that is the polarization developed in the z-direction per unit stress applied in the x-direction [19].

Model construction 2.2.1. Model
The FE model consisted of a simply supported plate with a piezoelectric patch attached on its surface. Six different patch shapes were modeled in order to study their performance concerning the change in the patch shape. The 1060 aluminum alloy plate had a thickness of 1 mm, while the piezoelectric patches had a thickness of 0.5 mm and all had the same surface area (i.e., same volume and mass). The piezoelectric patches were located at ¾ of both the plate's length and width, as shown in figure 3. This would allow capturing the bending of more bending modes within the specified frequency range. To induce randomness in the plate, 30 point masses were distributed randomly on 90% of the plate's inner area. The total mass of the point masses was 30% of the weight of the plate. The randomization of the point masses created uncertainties in the structure that induce the effect of infinite plate structures which in turn produced the response of Non-DS. Figure 3 presents the benchmark model adopted for this study. Figure 4 presents the six shapes that are considered in this study. The piezoelectric material of the piezoelectric patches used in the model is the ceramic type piezoelectric PZT-5A that is polarized along the z-axis (of figure 3). The material properties of the aluminum plate are compiled in table 1 below and those of the PZT-5A are compiled in table 2.

Element types and meshing
Four element types were used to construct the FE model: SOLID185 for the plate, SOLID5 for the piezoelectric patch, MASS21 for the point masses, and CIRCU94 for the external resistance. Information about the selected elements and their capabilities can be found in ANSYS ® theory reference [21]. SOLID185 is a 3D structural solid element and is defined by eight nodes having three translation degrees of freedom (DOFs) at each node. SOLID5 is a 3D coupled-field solid element that has structural and piezoelectric field capability with coupling between both fields. SOLID5 has eight nodes with up four DOFs at each node: three translation DOFs and electric potential DOF. The mesh of the FE model is illustrated in figure 5.

Boundary conditions
The FE model included a mechanical boundary condition on the plate and an electrical boundary condition on the piezoelectric patch. Since 3D elements only have 3 translation DOFs and do not have any rotational DOFs, modeling simply supported boundary conditions using these elements could be complicated. Therefore, the plate was simply supported by constraining the displacement DOFs of the lower edges of the plate, as shown in figure 3. This approach gives a similar boundary condition to the simply-supported boundary condition that is usually applied when considering plate elements. As for the electrical boundary condition, the nodes on the upper surface of the piezoelectric patch were coupled to have the same electric potential. Consequently, the surface nodes would act as one node, and the voltage would be the same throughout the surface, as shown in figure 3. Similarly, the nodes on the bottom surface of the patch were coupled to have the same electric potential; this is set to zero since this was considered as the electric ground. Figure 3 also shows the modeling of the electric circuit used in the FE model to measure the electric potential produced by the piezoelectric patch. The electric circuit was modeled by creating electrodes at the lower and upper surfaces of the patch. The nodes at the bottom surface were coupled and electrically grounded, whereas the nodes at the top surface were coupled to the same voltage. A high resistive load was then connected between the top and bottom electrodes of the piezoelectric patch, and ultimately, the electrical potential could be  obtained at the top electrode. Figure 6 shows a flowchart of the undertaken procedures performed throughout the FE modeling.

Results and analysis
MOF calculations of the proposed FE model were conducted to gain an insight into the most proper frequency range to use in the simulations. The MOF specifies the frequency at each mode which in return identifies the number of vibration modes within the selected frequency range. The Non-DS plate is typically subjected to structural uncertainties; therefore, one has to ensure that a proper frequency range is selected whereby the nondeterministic behavior can be observed. A 5% loss factor ( ) h was selected; this corresponds to a damping ratio ( ) z of 2.5%, since h z = 2 . The plate was subjected to a point force excitation of 10 N and the frequency range defined for the harmonic analysis was from 0-1200 Hz. Around 20 bending modes could be captured within the selected frequency range. The MOF was calculated using equations (1)-(3). The flexural rigidity was calculated using equation (3). A value of · 6.546 Pa m 3 of flexural rigidity was obtained and substituted in equation (2) to calculate the modal density. The modal density was retrieved from equation (2) was equal to / 0.0765 modes Hz and substituted in equation (1) to calculate MOF over the entire frequency range.
As for the natural frequency equation of the mn-th mode for simply-supported rectangular plates, it is expressed as: Where subscripts m and n refer to the half-wave number in the x and y directions. L x and L y are the length of the plate in the x and y directions. This equation was used to predict the amount of bending modes occurring within the specified frequency range. A modal analysis was performed prior to the harmonic analysis to identify and analyze the natural frequencies and mode shapes of the plate. Such analysis is crucial to visualize the mode shapes of the structure and ensure that the PEH will experience strain and thus effectively convert mechanical energy into electric power. For the harmonic analysis, 20 runs were performed for each piezoelectric patch shape. The 20 responses were then averaged to obtain the ensemble average response.

Wavelength calculations
In order to assure a maximum voltage generation, the plate's bending wavelength deformation was calculated and compared with the dimension of the six piezoelectric patches attached to the plate. The calculations were performed for the highest frequency (1200 Hz in this study) to obtain the shortest wavelength deformation. The plate's bending wavelength deformation was calculated using the following equation: Where D is the flexural rigidity, m p is the mass of the plate and w is the frequency in / rad s. The bending wavelength deformation of the plate was found to be equal to 8.63 cm and it was observed that the plate's wavelength deformation is above the size of the six piezoelectric patches.
Additionally, the number of elements per wavelength was calculated using: Where l m and l n are respectively the wavelengths in the x and y directions. At the maximum frequency of 1200 Hz (shortest wavelength), the mode shape function had the numbers of maximum = m 17 or = n 17. The wavelengths in the x and y directions were equal to 117.6 mmand 141.2 mm respectively. Thus, the commonly applied rule of thumb to use at least six elements per wavelength was practiced in this study. The values obtained from the equations above show that in order to have six elements per wavelength, the element size must be around 19.6 mm. In this study, the element size used to mesh the FE model was 12 mm.

Driving point response
To investigate the vibrational behavior of the structure, the driving point receptance was analyzed. Figure 7 shows the driving point receptance of the plate when attached to different piezoelectric patches (shapes and dimensions are in figure 4). For each piezoelectric patch, 20 simulations were conducted with random point masses attached to the plate. Figure 8 compiles the average driving point receptance of the six piezoelectric patches. Examining the driving point receptance clearly shows that the structure has low-variance in the lowfrequency range, whereas the variance is much higher in the mid-and high-frequency ranges (above 400 Hz). Figure 9 shows the electrical potential responses ( ) EP for the 20 randomized structures (for each of the six patch shapes). The peaks could be seen occurring at the natural frequencies of the structure, which were observed in figure 7. In the low-frequency range, the resonance peaks are distinct, whereas, in the high-frequency range, the resonance peaks are no longer distinct. At high frequencies, the length of the bending wavelength became comparable to the distances between the random point masses, and thus, the structure's response became more sensitive to structural uncertainties which in turn increases the variance. The response at the low-frequency range (i.e., low MOF) had high variance, while at the high-frequency range (high MOF), low response variance could be observed. The average response was smoothed out as the system vibrates at higher frequencies. The power dissipated over the plate's resistor was calculated for all the patches. The average power output was potted against the frequency and presented in figure 10.

Cumulative power response
Cumulative frequency analysis was performed to get more insight into the performance of the six patch shapes within different frequency ranges. Figure 11 depicts the cumulative frequency average power ( ) CFAP . It was obtained by calculating the integral of the power response ( ) P over the studied frequency range ( ) f as: From figures 11 and 12 it was observed that the quadrilateral patch exhibited the highest performance at low and high frequencies. The piezoelectric performance of the octagonal and triangular patches was 9% and 15% lower than that of the quadrilateral patch, respectively. It was also deduced that the performance of the octagonal patch was better than that of the triangular patch. However, at the frequency range 150-600 Hz, the triangular patch demonstrated a better efficiency, but it decreased as the system vibrates at higher frequency ranges. As for the dodecagonal and circular patches, they somewhat showed a similar performance, especially at a higher frequency. However, it was observed that the circular patch had the lowest overall piezoelectric performance at both low and high-frequency ranges. Additionally, the cumulative power average showed that more efficient piezoelectric performance could be  achieved as the system vibrated at a lower frequency. Overall, it could be deduced that the piezoelectric patches with an even number of edges possessed a better energy harvesting characteristic. Nevertheless, the performance of these patches declined by increasing the number of edges. The same case applies to patches with an odd number of edges.  The effect of the infinite plate structures can be further approached by inducing more randomness to the structure. The randomized condition of the plate structure creates the response uncertainties required to approach the effect of an infinite plate; that in turn, generates the response of a Non-DS. To promote further dynamic randomization, the structure was subjected to 20 time-harmonic randomly located point forces each of 5N. The 30 randomly located point masses are still attached to the plate. Figure 13 is an illustration of the randomly distributed point forces. Under this excitation condition, the performance of the proposed PEHs attached to the host plate was investigated by employing the ensemble average response. The ensemble responses were then analyzed to further investigate the effect of the piezoelectric patch shape on the developed voltage.     Figure 14 shows the ensemble electrical potential responses of the six proposed shapes. The average responses for the power output are plotted together in figure 15. Figure 16 shows the cumulative frequency average power CFAP of the six piezoelectric shapes. This plot gives more insight into the performance of the piezoelectric shapes at different frequencies. In the CFAP plot, it could be noticed that the performance of the piezoelectric patches was higher in the low-frequency range as the resonant peaks of the low-frequency range were more distinct and contributed to higher amplitudes. As for the overall performance, the dodecagonal patch exhibited the highest performance, followed by the triangular patch, with the latter contributing to higher performance within the low-frequency range. The circular patch had the lowest harvesting performance in the overall frequency range.

Conclusion
In this paper, the energy harvesting performance of piezoelectric patches with different shapes that are attached to a thin non-deterministic plate was investigated. Six shapes were investigated: triangular, quadrilateral, pentagonal,  octagonal, dodecagonal, and circular. The purpose of this investigation was to define the shape of the piezoelectric harvester that possesses the highest energy harvesting performance. A finite element model constructed in ANSYS ® was used to simulate the system at hand. The dynamics of the host plate were randomized by randomly distributing point masses on the plate to achieve the characteristics of the non-deterministic substructure. The plate was simply supported on all edges and subjected to a point force excitation. In addition, the modal overlap factor was analyzed in order to identify the proper frequency range to use. It must be noted that the orientation and location of the patches are important for low-frequency energy harvesting. However, for a high-frequency response of an uncertain structure, the orientation does not affect the ensemble average. In this work, the orientation was maintained for all patches where one side of the patch is aligned parallel with the plate's side. The electrical potential response developed across the piezoelectric patch was obtained and analyzed. It was found that the quadrilateral shape has the highest performance at various frequency ranges, while the circular patch is the lowest. However, an increase in the performance of the dodecagonal was noticed as the host plate was subjected to random point force distribution. Overall, it was deduced that the performance of piezoelectric harvesters could be optimized by using shapes with an even number of edges.

Data availability statement
The data that support the findings of this study are available upon reasonable request from the authors.