The Piezoelectric Phenomenon in Energy Harvesting Scenarios : A Theoretical Study of Viable Applications in Unbalanced Rotor Systems

The present paper deals with the promising energy harvesting applications of a composite piezoelectric metal support that is properly designed for the rotor of a mechanical system. The aim is to determine whether the vibrational power coming from the static residual imbalance, which is generally considered to be an undesired and useless side-effect of the rotation, can be converted into electric power and then stored to be used in other applications. The analysis, starting from the Jeffcott rotor model and the piezoelectric constitutive equations, has been carried out by developing an approximated linear model of a piezoelectric support, in order to theoretically evaluate the performance and the feasibility of the proposed system. The accuracy of the exploited analytical model has been validated for both static and dynamic operations by 3D Ansys® Mechanical APDL. Finally, a MatLab®/Simulink® model has been built to simulate the electric behavior of the piezoelectric material, and to estimate the power that it is possible to extract via an alternative/direct current converter (AC/DC converter). The numerical results achieved confirm the effectiveness of the proposed energy-harvesting system.


Introduction
Energy harvesting based on piezoelectric technology has been developed over the last decades for charging portable electronic devices and/or for micro-sensors and actuators, where the prominent aspect is to ensure the operation of the electronic devices [1,2].Indeed, a literature survey on piezoelectric energy harvesters shows that they are classified into three groups with reference to their correspondent sizes: macro-and mesoscale, MEMS (Micro Electro-Mechanical Systems) scale, and nanoscale [3].However, this approach seems that neglect the possibility for employing this technology in order to scavenge and store energy [4].In an actual energetic context where the global demand is always increasing [5,6], it is needful to exploit any source of renewable energy, in order to improve eco-sustainability.This work is focused on the extraction principle of a form of energy that is usually available in any rotating system: side vibrations.Indeed, via the piezoelectric effect, which is the capability of some solid materials to convert mechanical energy into electric energy and vice versa by accumulating electric charge in response to mechanical stress, it is possible to exploit these vibrations, and, for instance, deliver the extracted energy to a storage device.
Piezoelectric materials are historically associated with micro-actuators/sensors, or, in the last two decades, micro-power generation [7][8][9][10][11].The main reasons that have led to their establishment in these sectors are the simplicity of their shaping, especially for piezocomposites such as PZT (Lead zirconate titanate), and the impressive performances that are possible to be obtain with their very small dimensions [12][13][14][15][16][17].In this context, the accuracy with which it is possible to predict the mechanical/electric responses of these materials has been the key factor to their success.In the last decade, in particular, the attention of the scientific community has been mainly focused on applications other than in the sensor field, i.e., the energy harvesting applications.
This paper aims to demonstrate that the piezoelectric phenomenon can be exploited, without performance losses, in order to scavenge the energy related to the mechanical vibrations of a system where the main power is delivered through rotation.

Materials and Methods
The proposed application employs the direct piezoelectric effect [18,19]: the piezoelectric material is driven by mechanical vibrations coming from an unbalanced rotor, and the consequently produced voltage is collected through an AC/DC converter, so that when a load is connected, electrical power can be generated and/or stored.The motor concept can be seen in Figure 1a.The piezoelectric material is inserted in a piezoelectric support, modelled in a unimorph-like fashion [12], leading to a double-layered piezoelectric-steel beam.In particular, each of the four positioned supports is embedded at one end and rigidly connected to the ball bearings at the other end, as schematized in Figure 1b.Naturally, in contrast to the motor-operating conditions [18], no additional electric sources are required.
Energies 2018, 11, x FOR PEER REVIEW 2 of 19 mechanical/electric responses of these materials has been the key factor to their success.In the last decade, in particular, the attention of the scientific community has been mainly focused on applications other than in the sensor field, i.e., the energy harvesting applications.This paper aims to demonstrate that the piezoelectric phenomenon can be exploited, without performance losses, in order to scavenge the energy related to the mechanical vibrations of a system where the main power is delivered through rotation.

Materials and Methods
The proposed application employs the direct piezoelectric effect [18,19]: the piezoelectric material is driven by mechanical vibrations coming from an unbalanced rotor, and the consequently produced voltage is collected through an AC/DC converter, so that when a load is connected, electrical power can be generated and/or stored.The motor concept can be seen in Figure 1a.The piezoelectric material is inserted in a piezoelectric support, modelled in a unimorph-like fashion [12], leading to a double-layered piezoelectric-steel beam.In particular, each of the four positioned supports is embedded at one end and rigidly connected to the ball bearings at the other end, as schematized in Figure 1b.Naturally, in contrast to the motor-operating conditions [18], no additional electric sources are required.PZT, grey layer: steel).

Piezoelectric Material Characterization
The employed piezoelectric material for this application is PZT, a lead-titanium-zirconium alloy belonging to the symmetry class  21 [20][21][22][23][24][25][26].In particular, PZT 850, produced by the American group Piezo Ltd., has been considered.This material, whose datasheet [26] values are summed up in Table 1, belongs to the Navy Type II [27], and has been chosen due to its high piezoelectric coefficients d 3i and its good deformability.

Piezoelectric Material Characterization
The employed piezoelectric material for this application is PZT, a lead-titanium-zirconium alloy belonging to the symmetry class C 2v 21 [20][21][22][23][24][25][26].In particular, PZT 850, produced by the American group Piezo Ltd., has been considered.This material, whose datasheet [26] values are summed up in Table 1, belongs to the Navy Type II [27], and has been chosen due to its high piezoelectric coefficients d 3i and its good deformability.The following form of the constitutive equations, among the four possible forms [23], has been employed: The relations (1), derived from the differentiation of the internal energy of a piezoelectric crystal, require the definition of the tensors s E ij , e in and ε T nm , while the other mechanical, electric, and piezoelectric tensors can be consequently derived by using the relations contained in Table 2.For single-term characterization, the relations contained in [28,29] have been employed.

Elastic Coefficients
Piezoelectric Coefficients Dielectric Coefficients By using an orthotropic model for the anisotropic PZT material, and by putting K T 11 = 1980 and ν 31 = 0.39, the compliance tensors s E ij and e in are derived:

Analytical Model
Since the four piezoelectric devices are arranged symmetrically in two couples of 90 • -shifted supports, the system behavior can be derived by characterizing the mechanical stress of a generic couple of supports.Table 3 particularizes the mechanical stress (force or torque) for each support (lower and upper), depending on the phase shift θ of the rotor mass center with respect to the y axis.Indeed, for each θ, two kinds of mechanical stress can be identified: a buckling/traction stress (driven by the inertia force component N along the support axis) and a bending one (driven by the torque T, associated with the complementary inertia force component).The double-layered support has been modeled for static and dynamic analyses as a one-dimensional beam [30,31] with only one degree of freedom, i.e., the displacement of the free end in the direction orthogonal to the axis of the beam.The analysis was performed with reference to one of the lower supports.Naturally, the relative results can be easily extended to the upper supports.Moreover, it is assumed that the rotor-generated mechanical stress is equally distributed and concentrated on the supports' free ends, whose masses are neglected.In this context, the PZT material is assumed to be isotropic and characterized by a Young constant: E PZT = 74 GPa.The analytical analysis was validated by a numerical one in the Ansys ® environment, where the actual orthotropic model of the material and the relative constitutive equations were implemented.The geometric and mechanical characteristics of the beam are listed in Table 4, while the mono-dimensional scheme is shown in Figure 2. The double-layered support has been modeled for static and dynamic analyses as a onedimensional beam [30,31] with only one degree of freedom, i.e., the displacement of the free end in the direction orthogonal to the axis of the beam.The analysis was performed with reference to one of the lower supports.Naturally, the relative results can be easily extended to the upper supports.Moreover, it is assumed that the rotor-generated mechanical stress is equally distributed and concentrated on the supports' free ends, whose masses are neglected.In this context, the PZT material is assumed to be isotropic and characterized by a Young constant: E PZT = 74 GPa.The analytical analysis was validated by a numerical one in the Ansys ® environment, where the actual orthotropic model of the material and the relative constitutive equations were implemented.The geometric and mechanical characteristics of the beam are listed in Table 4, while the monodimensional scheme is shown in Figure 2.  The amplitude of the force exerted on the free end has been calculated from the bending stress, constrained to its maximum admissible value T MAX .By imposing the ratio T ADM T MAX to be equal to 1.25, the target stress force is found: The mechanical stress due to the normal loads (compressive/tensile), which represents 0.08% of the total axial stress, has been neglected.
While in the static analysis the force (4) generates a constant displacement along the z-axis of the beam, in the dynamic one, it represents the amplitude of an alternative stress that produces a motion of the free end of the beam.In each case, the relation between the free-end displacement and the PTZ voltage will be derived.

Static Analysis: Free-End Displacement
With reference to Figure 2, the free-end displacement is computed by the Virtual Work Principle (VWP).In particular, the strain field S y,z and the displacement field v(z,y) along the beam axis are linked by the relation: S y,z = ∂v z,y ∂y In a generic system Π, the displacement field can be determined by minimizing the relative potential energy with the Rayleigh-Ritz method: the displacement field is therefore a linear The amplitude of the force exerted on the free end has been calculated from the bending stress, constrained to its maximum admissible value T MAX .By imposing the ratio T ADM T MAX to be equal to 1.25, the target stress force is found: The mechanical stress due to the normal loads (compressive/tensile), which represents 0.08% of the total axial stress, has been neglected.
While in the static analysis the force (4) generates a constant displacement along the z-axis of the beam, in the dynamic one, it represents the amplitude of an alternative stress that produces a motion of the free end of the beam.In each case, the relation between the free-end displacement and the PTZ voltage will be derived.

Static Analysis: Free-End Displacement
With reference to Figure 2, the free-end displacement is computed by the Virtual Work Principle (VWP).In particular, the strain field S(y, z) and the displacement field v(z, y) along the beam axis are linked by the relation: In a generic system Π, the displacement field can be determined by minimizing the relative potential energy with the Rayleigh-Ritz method: the displacement field is therefore a linear combination of a function f(y, z), and its m − 1 differentials via m coefficients are calculated as follows: Energies 2019, 12, 708 5 of 20 In the case of the bending beam, according to the Euler-Bernoulli theory, the axial displacement field can be given as a function of the rotation ϕ x and the orthogonal displacement field w(y): v(y, z) = zϕ x = z ∂w(y) ∂y (7) with w(y) being a Hermite polynomial function.Since the axial strain field is obtained as the differential of the axial displacement field v(y, z) with respect to y, the z displacement field has to be a second-degree polynomial, so that the strain field is at least constant.However, when a continuous domain is discretized via finite elements that are linked through several nodes, a displacement field is valid if there is a correspondence between the total degree of freedom of the nodes and the number of m coefficients used.Therefore, the minimum degree possible for the sought polynomial is the third degree: Given that y ∈ [0, l] and the boundary conditions w(0)= 0 and .w(0)= 0 imposed by the joint configuration, it can be simply verified that the first two constants are null, so that the elastic internal energy is: Naturally, E el is related to the displacement of the beam free end by the relation: From ( 9) and ( 10), c 2 and c 3 are determined.In particular: Once the shape function is introduced: the strain field can be written as: Concerning the (12), which refers to a homogeneous beam, an equivalent expression can be given for a composite beam, using an approximated calculation of the stiffness [27]: The modeled beam is a composite, so that ( 15) represents an approximation with respect to the Classical Lamination Theory.The comparison between the results obtained by the analytical model and the finite element (FE) model will give a measure of this approximation.The parameter D is the Energies 2019, 12, 708 6 of 20 equivalent to the Young constant of a homogenous beam multiplied by its bending inertial moment, so that for a composite beam, it results in the following: The inertial moment of the generic layer is related to its thickness, its width, and the orthogonal distance between its center of mass and the neutral bending axis.The neutral bending axis position can be found according to the geometric/mechanical properties of the layers considered.Assuming the lower layer as the steel one, the neutral bending axis coordinate is: Consequently, it results in the following: An equivalent stiffness of the PZT material can be determined by using the constitutive piezoelectric Equation (1).Since: by assuming a mono-dimensional and constant electric field, and by imposing V A = 0, the voltage V B is obtained: The electric stiffness of the system can be calculated by imposing that the internal work L int is the result of the work where the stresses perform for the virtual strains δS, and the one that the electric displacement performs for the virtual electric field δE.Thus, by expressing the stress and the electric displacement according to Equation (1), the following relation can be written: Using ( 14) and ( 23), the integrals in ( 24) can be evaluated by assuming a linear behavior, and by splitting each contribution in the correspondent layer.In particular, the steel layer will contribute only to the stiffness term, while the PZT one will contribute to all terms.Given the isotropic assumption for the PZT and that the Young modules are constants, the mechanical stiffness is calculated as: where the term g(y) = g(z, y)/z is defined in order to separate the surface and the linear integrals.The surface integral of z 2 is the inertial moment with respect to the x-axis, and it has already been calculated in (18).By computing the remaining integral, the mechanical stiffness of the beam is derived: which obviously coincides with (20), since this term is not affected by the electro-mechanical effect.
Energies 2019, 12, 708 7 of 20 The volume piezoelectric coupling constant is given by: Based on the computed strain along y and electric field along z, by taking into account the piezoelectric constant e 32 = −10.1 NV/m, (27) becomes: Hence, the PZT volume dielectric constant is calculated as: The polarization occurs along z, and it is assumed to be independent of the length and the width of the beam.The remaining integral can be therefore computed as: Analogously, the external work is based on the contributions of the work surface charge density σ e , and the external force F: By constraining L ext to zero, the integral constitutive equation for the beam is derived: from which the electric equivalent stiffness of the beam is obtained as: Equation ( 33) is higher than the (20).Indeed, the electric properties of PZT tighten up the system and, consequently, the actual displacement of the free end is lower than (21): The deformed shapes of the beam have been represented in Figure 3, solving the displacement field in a Matlab ® script.
Equation ( 33) is higher than the (20).Indeed, the electric properties of PZT tighten up the system and, consequently, the actual displacement of the free end is lower than (21): The deformed shapes of the beam have been represented in Figure 3, solving the displacement field in a Matlab ® script.After the analytical solution, a 3D model of the support has been produced with Ansys ® Mechanical APDL, to determine both the mechanical and the electric-mechanical responses, in terms of the free-end displacement of the beam.SHELL281 elements have been used in purely mechanical analyses, with these being the most suitable for the analysis of a thin composite structure.The deformed shapes for the isotropic and the orthotropic models are depicted, respectively, in Figure 4.After the analytical solution, a 3D model of the support has been produced with Ansys ® Mechanical APDL, to determine both the mechanical and the electric-mechanical responses, in terms of the free-end displacement of the beam.SHELL281 elements have been used in purely mechanical analyses, with these being the most suitable for the analysis of a thin composite structure.The deformed shapes for the isotropic and the orthotropic models are depicted, respectively, in Figure 4.After the analytical solution, a 3D model of the support has been produced with Ansys ® Mechanical APDL, to determine both the mechanical and the electric-mechanical responses, in terms of the free-end displacement of the beam.SHELL281 elements have been used in purely mechanical analyses, with these being the most suitable for the analysis of a thin composite structure.The deformed shapes for the isotropic and the orthotropic models are depicted, respectively, in Figure 4. Comparing the analytical results (21) with the simulation ones in 4a and 4b, and (34) with 4c and 4b, the deviations in Table 5 have been found: With reference to the first column of the table, it is evident that in an isotropic condition, the FE model is stiffer than the analytical one.This aspect can be attributed to the discretization operated by the meshing process.Naturally, the orthotropic model registers a higher deviation, even though the accuracy is still acceptable.A purely mechanical analysis has been performed with the SOLID226 model, whose results are consistent with analytical ones.
All of the results from the analytical analyses are summarized in Table 6.It should be pointed out that the Young constant derived from the piezoelectric results of the SOLID226 model gives E PZT = 72.7 GPa.

Static Analysis: Buckling
When a structure has a preponderant dimension with respect to the other two, an instability may occurs if the mechanical stress acts along the considered direction.To analysis this aspect in the case of an embedded beam with an axial force exerted on the free end, the beam equilibrium is considered with respect to the deformed shape, since the hypothesis of small displacements no longer applies.
The critical buckling load for a composite beam can be evaluated by the Euler relation ( 16): where l is the beam length and C is a constrain factor.The evaluated critical load is two orders of magnitude bigger than the load considered in this context, F = 54.88N, so that the stability test is satisfied.By considering the slim ratio l/k, the buckling analysis also has to be extended to the Euler buckling hyperbole, which is defined by the relation: Finally, the Johnsons critical stress also has to be considered: The intersection between (36) and the (37), represented in Figure 5a, gives the threshold slim ratio, such that if the slim ratio of the considered beam is smaller than this limit, the critical stress has to be calculated with (36), otherwise (37) must be used.In this case, the threshold slim ratio is higher than the beam slim ratio, so that the critical load is evaluated via (37): From (38), it can be deduced that the bucking verification is also satisfied.The same analysis has been performed on the SHELL281 model, and the eigenvalues in Figure 5b have been obtained.The buckling analysis can indeed be expressed as an eigenvalues problem resulting from the deformed shape equilibrium: where K is the stiffness matrix of the system, K G is the geometric stiffness matrix of the system, whose value depends on the baseline load fixed P * , and λ is an eigenvalue such that the characteristic polynomial |K + λK G | = 0. Supposing a baseline load P * , a static analysis followed by the buckling analysis has been performed, obtaining the eigenvalues in Figure 5b.It has to be noted that since the beam is a composite, the instability could be misinterpreted as an excessive bending reaction.With respect to the first positive eigenvalue, the Euler critical load can be determined as: P cr,mod = 49.8P b = 2733 N (40) so that the deviation between the analytical model ( 35) and the FE model ( 40) can be calculated as: The software has therefore overestimated the structural resistance, although the deviation is acceptable.38), it can be deduced that the bucking verification is also satisfied.The same analysis has been performed on the SHELL281 model, and the eigenvalues in Figure 5(b) have been obtained.The buckling analysis can indeed be expressed as an eigenvalues problem resulting from the deformed shape equilibrium: (K+λK G )v = λP * (39) where K is the stiffness matrix of the system, K G is the geometric stiffness matrix of the system, whose value depends on the baseline load fixed P * , and λ is an eigenvalue such that the characteristic polynomial |K+λK G | = 0.
Supposing a baseline load P*, a static analysis followed by the buckling analysis has been performed, obtaining the eigenvalues in Figure 5(b).It has to be noted that since the beam is a composite, the instability could be misinterpreted as an excessive bending reaction.With respect to the first positive eigenvalue, the Euler critical load can be determined as: P cr,mod = 49.8P b =2733 N (40) so that the deviation between the analytical model ( 35) and the FE model (40) can be calculated as: The software has therefore overestimated the structural resistance, although the deviation is acceptable.

Dynamic Analysis: Fatigue
The mechanical stress applied on the free end of the support is an alternative bending with amplitude: The frequency of the load is equal to the speed of the rotor, set to 3000 rpm, 4500 rpm, and 6000 rpm.The amplitude is maintained as a constant for the variation of the speed, by adjusting the eccentricity value.In this case, the term eccentricity has been improperly used to represent the product between the mass and the actual eccentricity, i.e., the distance between the center of mass and the center of rotation.The three load curves have been plotted in Figure 6a with respect to the period of the slower one, while the Whӧler plot of the material has been represented in Figure 6b.This plot is divided in three successive parts: the first part, low fatigue cycle (LFC), starts from the ultimate strength of the material, and the third part is a plateau representing the stress magnitude that corresponds to an infinite life of the component (Cutoff), while for the second part high fatigue cycle (HFC) the following relation has been used: ln S f = ln a+b ln N (43)

Dynamic Analysis: Fatigue
The mechanical stress applied on the free end of the support is an alternative bending with T a = T max −T min 2 = 0.24 10 8 GPa (42) The frequency of the load is equal to the speed of the rotor, set to 3000 rpm, 4500 rpm, and 6000 rpm.The amplitude is maintained as a constant for the variation of the speed, by adjusting the eccentricity value.In this case, the term eccentricity has been improperly used to represent the product between the mass and the actual eccentricity, i.e., the distance between the center of mass and the center of rotation.The three load curves have been plotted in Figure 6a with respect to the period of the slower one, while the Whöler plot of the material has been represented in Figure 6b.This plot is divided in three successive parts: the first part, low fatigue cycle (LFC), starts from the ultimate strength of the material, and the third part is a plateau representing the stress magnitude that corresponds to an infinite life of the component (Cutoff), while for the second part high fatigue cycle (HFC) the following relation has been used: with S f(10 3 ) as the maximum stress amplitude that can be applied for 10 3 cycles before collapse, and S e as the cutoff stress limit.While S f(10 3 ) can be estimated from the ultimate strength via a factor f < 1, which decreases as the ultimate strength increases, S e can be assigned once the stress amplitude S e is identified.Once S e , the fatigue limit of a sample that is stressed with alternative bending, is estimated as being half the ultimate strength of the material, S e is obtained as: The Whöler plot parameters and the service life of the component can then be evaluated: Inter alia, in this case, since the cutoff limit ( 45) is above the stress amplitude exerted on the material (42), the service life of the component has been assumed as being infinite.where a and b are two parameters derived by forcing the line (43) to pass through the following points: A (10 3 ,S f (10 3 )) B (10 6 ,S e ) (44) with S f(10 3 ) as the maximum stress amplitude that can be applied for 10 3 cycles before collapse, and S e as the cutoff stress limit.While S f(10 3 ) can be estimated from the ultimate strength via a factor f<1, which decreases as the ultimate strength increases, S e can be assigned once the stress amplitude S e ' is identified.Once S e ', the fatigue limit of a sample that is stressed with alternative bending, is estimated as being half the ultimate strength of the material, S e is obtained as: S e a k b k c d k e S e ' =241 MPa (45) The Whӧler plot parameters and the service life of the component can then be evaluated: Inter alia, in this case, since the cutoff limit ( 45) is above the stress amplitude exerted on the material (42), the service life of the component has been assumed as being infinite.

Dynamic Analysis: Equation of Motion
The rotor-support system can be modeled according to the Jeffcott rotor model, assuming a rigid rotor suspended on elastic supports.The rotor is assumed to be fixed at the mid-section of the drive, so that the system is symmetrical with respect to the yz plane.The four supports, assumed to be identical, have been modeled as a spring-damper system with uniform stiffness and dumping coefficients in the y and z directions.If the rotor is statically unbalanced, the motion equations can be written as [28]: with m as the rotor mass and e as the eccentricity.Given the aforementioned assumptions, the motions along the two axis are therefore decoupled.In particular, the force amplitude evaluated in ( 4) is exerted on each support, so that the following relation applies: 2F = meω 2 , (48) Given that (48) is constant, the three resulting eccentricity values are immediately obtained; they are listed in Table 7.

Dynamic Analysis: Equation of Motion
The rotor-support system can be modeled according to the Jeffcott rotor model, assuming a rigid rotor suspended on elastic supports.The rotor is assumed to be fixed at the mid-section of the drive, so that the system is symmetrical with respect to the yz plane.The four supports, assumed to be identical, have been modeled as a spring-damper system with uniform stiffness and dumping coefficients in the y and z directions.If the rotor is statically unbalanced, the motion equations can be written as [28]: m ..
with m as the rotor mass and e as the eccentricity.Given the aforementioned assumptions, the motions along the two axis are therefore decoupled.In particular, the force amplitude evaluated in ( 4) is exerted on each support, so that the following relation applies: Given that (48) is constant, the three resulting eccentricity values are immediately obtained; they are listed in Table 7.Clearly, the eccentricity decreases with progressively smaller decrements as the speed increases.By considering a rotor with a length of 80 mm and a diameter 50 mm, a balanced mass, m 1 = 1.78 kg, is derived.The eccentricity property is then conferred to the system removing n identical masses with a constant phase shift of 15 • , starting from the mass m n at θ = 0 • .Since the masses are symmetrically positioned with respect to the y-axis, the resultant mass center is simply shifted along the y-axis.By putting n = 13, and by solving the corresponding system of equations, the values listed in Table 8 are calculated: The rotor masses included in Table 8 correspond to twice the mass of the rotor supports, which, indeed, has been schematized as beams stressed by a concentrated force applied at the mid-section.This approximation leads to an underestimation of the resonance frequencies of the system.
With respect to the system stiffness, each support contributes to the values of the analytic electric stiffness and the SOLID226 electric stiffness found in Table 6.
Assuming the stiffness constant with the speed, both the resonance frequencies and the critical dumping coefficients can be evaluated according to the following relations (see also Table 9): As expected, as the mass increases at a constant stiffness, and natural frequency decreases.
Energies 2019, 12, 708 13 of 20 The damping action of the system represents the hysteretic/friction phenomena that occur within and between the system parts.To calculate the corresponding coefficients, various methods can be used, and generally, the experimental interpolation of the data collected from the actual system methods are preferred.In this case, however, the proportional analytical damping model has been adopted: Once the coefficients are assumed to be equal, the corresponding value can calculated by fixing the magnitude of the damping ratio ξ.Inter alia, if ξ < 1 is fixed, the system law of motion will be an actual vibration, otherwise it would return with an exponential behavior to its initial condition.This assumes that the system has not been irreparably perturbed, in which case, it can either reach a new equilibrium state, or remain in a non-equilibrium state that can eventually take it to collapse.The values of the damping coefficients are reported in Table 10 for different system configurations.From the parameter characterization, it is finally possible to determine the motion law.By focusing on an upper support (the same results can be extended to a lower support), the dynamic Equation (47) can be written as: from which: The value of the gain factor χ, of the speed ratio β = ω ω n and of the phase angle ψ are reported in Table 11 for different speed values.As expected, the gain factor is maximum when the rotor speed equals 3000 rpm, as the resonance frequency of the system is at 3190 rpm.It should be noted that the piezoelectric materials guarantee the best electro-mechanical conversion at around their resonance frequency, where indeed, the impedance of the correspondent equivalent electrical flow is minimized.The time behavior of motion, speed and acceleration at the three considered speeds is plotted in Figure 7.
1.48 1.49 -1.09 -1.08 -0.78 -0.77 As expected, the gain factor is maximum when the rotor speed equals 3000 rpm, as the resonance frequency of the system is at 3190 rpm.It should be noted that the piezoelectric materials guarantee the best electro-mechanical conversion at around their resonance frequency, where indeed, the impedance of the correspondent equivalent electrical flow is minimized.The time behavior of motion, speed and acceleration at the three considered speeds is plotted in Figure 7.The dynamic analyses have been performed on the SHELL281 FE model by including a MASS21 element whose RCs has been derived from the mass values shown in Table 9.With respect to the values listed in Table 8, the following deviations regarding the natural frequencies have been calculated: It should be pointed out that all the estimated errors are very small.

Electric Analysis
The charge Q that accumulates on the PZT surfaces as an effect of the applied mechanical stress can be collected by connecting a load, such as a resistor R L , via two electrodes anchored on the PZT surfaces.The current in the PZT-R L circuit can naturally be expressed as: By discretizing (53) with respect to a sample time interval ∆t, the following relation is obtained: Equation ( 54) can be solved once the initial value Q 0 is determined.Naturally, the evaluation of the Q n+1 is an iterative process, whose number of iterations depends on the ratio ∆t/R L .The voltage V can be computed, starting from the expression of the electric displacement flux density through the piezoelectric material: with q c as the surface charge.The electric displacement can be also expressed as a combination of the strain and the electric field according to the constitutive Equation (1), where, since only the component of the strain along the support axis applies, can be characterized by the piezoelectric coefficient e 3,2 : By imposing Q 0 = 0, by taking into account (14), and by expressing the free-end displacement according to (52), (56) leads to: The voltage on the upper surface of the PZT layer is therefore linear along the length of the beam, being maximum at y = 0, where the strain also is maximum.From (57), it can be deduced that when ∆t/R L ≤ 10, the voltage value is not affected by variations in ∆t or R L , while when ∆t/R L > 10, as the resistance value increases, so does do voltage values up to ∆t/R L ≥ 10 4 , at which point the voltage switches its sign and the solution becomes divergent.The aforementioned considerations are endorsed by Figure 8a,b, where the voltage is plotted in correspondence with y = l/2.The voltage behavior along the support axis is depicted instead in Figure 8c.

Equivalent Electric Circuit of PZT
From an electrical point of view, the mechanical behavior of a piezoelectric material can be modeled as an RLC series circuit [32,33]:

Equivalent Electric Circuit of PZT
an electrical point of view, the mechanical behavior of a piezoelectric material can be modeled as an RLC series circuit [32,33]: where inductance L M is half the rotor mass of the rotor, and resistance R M relates to the damping mechanical coefficient, while the electric stiffness corresponds to the inverse of the capacitance C M .Analogously, the piezoelectric behavior can be simply modeled by a capacitor C P series connected to a current generator that is controlled by the output of the aforementioned equivalent electric circuit.C P can be calculated as: The overall circuit is shown in Figure 9a, while all of the parameters are listed in Table 12.The parameters α and β result from the resonance frequency of the open-circuit configuration f 0 , which accounts for an equivalent capacitance that is given by the series C By putting the open-circuit resonance frequency at a slightly higher value than the short-circuit resonance frequency, the sought product of the parameters is obtained, so that by imposing α/β = 50 [30], their values are derived.
From Figure 9b, it can be noted that the voltage output of the piezoelectric layer is not perfectly sinusoidal, but more square-wave like.The voltage peak is around 153.5 V.  From Figure 9b, it can be noted that the voltage output of the piezoelectric layer is not perfectly sinusoidal, but more square-wave like.The voltage peak is around 153.5 V.The PZT surfaces are connected to an AC/DC converter, so that a continuous power P OUT can be produced and then stored by exploiting the piezoelectric effect.The designed converter is a full-wave single-phase rectifier.An external capacitance C C has been added in order to regulate the conduction mode, as represented in Figure 10a.The time behavior of the rectifier output voltage V L (t) is shown in Figure 10b, while that of the correspondent load current is depicted in Figure 10c.The extracted active power (see Figure 10d Given that the mechanical input power related to the vibration phenomenon can be calculated as:

System of Collection and Conversion of the Charge
the piezoelectric conversion efficiency can be finally derived: η PIEZO = P el P mecc =0.029 (64)

Results
The aforementioned analysis has been carried out for a single support, and the results obtained If the value of the capacitance C C is low, the converter works in continued conduction mode, and the average voltage value is given by the relation: If, on the other hand, the value of the capacitance C C is sufficiently high, the rectifier works in discontinuous conduction mode, and the output voltage can be approximated with the input voltage peak.
The time behavior of the rectifier output voltage V L (t) is shown in Figure 10b, while that of the correspondent load current is depicted in Figure 10c.The extracted active power (see Figure 10d) is equal to: P OUT = I av V av = 0.11 W (62) Given that the mechanical input power related to the vibration phenomenon can be calculated as: the piezoelectric conversion efficiency can be finally derived: η PIEZO = P el P mecc = 0.029 (64)

Results
The aforementioned analysis has been carried out for a single support, and the results obtained are summarized in Table 13: Given that the application employs four series-connected supports, the results of the whole system are listed in Table Table 14.Summary of the results concerning the entire application.

Conclusions
In this paper, an energy-harvesting system based on the piezoelectric effect has been proposed.The discussed system, based on properly arranged 90 • -shifted supports built upon a PZT layer and a steel one, is able to extract the vibration energy that is associated with the revolution of a statically unbalanced rotor.
The harvesting device has been schematized by a mono-dimensional isotropic analytical model, which has been validated in the Ansys ® Mechanical APDL environment and double-checked with Matlab ® .Furthermore, the time behavior of the main electric quantities has been studied by using the Simulink ® Matlab ® tool.
In particular, the proposed system is able to produce electric power of around 0.44 W with 3% efficiency.Similar applications, which employ the piezoelectric effect to produce mechanical power, are generally characterized by better conversion ratios.Indeed, the proposed implementation only accounts for coupling between axial stress and thickness polarization, consequently

Figure 1 .
Figure 1.(a) Ferroelectric rotor-piezoelectric supports (yellow) system with ball bearings (black) on the fixed base, as realized in SolidWorks; (b) Piezoelectric support scheme with the local coordinate system (yellow layer:

Figure 1 .
Figure 1.(a) Ferroelectric rotor-piezoelectric supports (yellow) system with ball bearings (black) on the fixed base, as realized in SolidWorks; (b) Piezoelectric support scheme with the local coordinate system (yellow layer: PZT, grey layer: steel).

Figure 2 .
Figure 2. Mono-dimensional scheme of the support with the local coordinate system and the mass corresponding to half the rotor mass concentrated at the free end.

Figure 2 .
Figure 2. Mono-dimensional scheme of the support with the local coordinate system and the mass corresponding to half the rotor mass concentrated at the free end.

Figure 3 .
Figure 3. (a) Deformed shape of the PZT-Steel composite beam; (b) Close up view.

Figure 3 .
Figure 3. (a) Deformed shape of the PZT-Steel composite beam; (b) Close up view.

Figure 3 .
Figure 3. (a) Deformed shape of the PZT-Steel composite beam; (b) Close up view.

Figure 4 .
Figure 4. Displacement vector of the composite beam for the isotropic mechanical model (a), isotropic electro-mechanical model (b), orthotropic mechanical model (c), orthotropic electro-mechanical model (d).Units are in mm, kg, N, Vmm N , MPa, Hz.
ln S f = ln a + b ln N (43) Energies 2019, 12, 708 11 of 20where a and b are two parameters derived by forcing the line (43) to pass through the following points:

Figure 7 .
Figure 7. Laws of motion, speed, and acceleration of the free end of the upper support at 50 Hz (a), at 75 Hz (b), and at 100 Hz (c).Frequency response of the beam according to Ansys ® evaluations (d).

Figure 8 .
Figure 8. Voltage-time behavior on the upper surface of the PZT layer at different resistance values in correspondence with ∆t/R L >10 (a) and ∆t/R L <10 (b).Voltage-space behavior along the support axis (c).Voltage-time behavior at 50 Hz, 75 Hz, and 100 Hz on the upper surface of the PZT layer (d).

Figure 8 .
Figure 8. Voltage-time behavior on the upper surface of the PZT layer at different resistance values in correspondence with ∆t/R L > 10 (a) and ∆t/R L < 10 (b).Voltage-space behavior along the support axis (c).Voltage-time behavior at 50 Hz, 75 Hz, and 100 Hz on the upper surface of the PZT layer (d).

Figure 9 .
Figure 9. Simulink ® model for the PZT equivalent electric circuit (a) and voltage output of the PZT layer (b).

Figure 9 .
Figure 9. Simulink ® model for the PZT equivalent electric circuit (a) and voltage output of the PZT layer (b).

Table 2 .
Relations between the electric, piezoelectric, and mechanic tensors.

Table 3 .
Main mechanical stress for the two supports.

Table 4 .
Composite support-beam geometric and mechanic properties.

Table 4 .
Composite support-beam geometric and mechanic properties.

Table 6 .
Static analysis results.

Table 7 .
Static analysis results.

Table 7 .
Static analysis results.

Table 8 .
Mass and eccentricity values for each rotor speed.

Table 9 .
Mass and eccentricity values for each rotor speed.

Table 10 .
Damping coefficients and damping rotations for the three speeds considered for the analytical model and the FE model.

Table 11 .
Amplification factors and phase angles of the law of motions for the three considered speeds.

Table 12 .
Equivalent electric circuit parameters.

Table 13 .
Summary of the results concerning a single support.