Performance Correction and Parameters Identification Considering Non-Uniform Electric Field in Cantilevered Piezoelectric Energy Harvesters

In the current electromechanical model of cantilevered piezoelectric energy harvesters, the assumption of uniform electric field strength within the piezoelectric layer is commonly made. This uniform electric field assumption seems reasonable since the piezoelectric layer looks like a parallel-plate capacitor. However, for a piezoelectric bender, the strain distribution along the thickness direction is not uniform, which means the internal electric field generated by the spontaneous polarization cannot be uniform. In the present study, a non-uniform electric field in the piezoelectric layer is resolved using electrostatic equilibrium equations. Based on these, the traditional distributed parameter electromechanical model is corrected and simplified to a practical single mode one. Compared with a traditional model adopting a uniform electric field, the bending stiffness term involved in the electromechanical governing equations is explicitly corrected. Through comparisons of predicted power output with two-dimensional finite element analysis, the results show that the present model can better predict the power output performance compared with the traditional model. It is found that the relative corrections to traditional model have nothing to do with the absolute dimensions of the harvesters, but only relate to three dimensionless parameters, i.e., the ratio of the elastic layer’s to the piezoelectric layer’s thickness; the ratio of the elastic modulus of the elastic layer to the piezoelectric layer; and the piezoelectric materials’ electromechanical coupling coefficient squared, k312. It is also found that the upper-limit relative corrections are only related to k312, i.e., the higher k312 is, the larger the upper-limit relative corrections will be. For a PZT-5 unimorph harvester, the relative corrections of bending stiffness and corresponding resonant frequency are up to 17.8% and 8.5%, respectively. An inverse problem to identify the material parameters based on experimentally obtained power output performance is also investigated. The results show that the accuracy of material parameters identification is improved when considering a non-uniform electric field.


Introduction
Energy harvesting is indeed the process of capturing and converting ambient sources of energy, such as light, heat, vibration, or motion, into usable electrical energy.It is a promising technology used to develop self-sustaining systems so that the need for battery replacement and disposal can be minimized [1][2][3].In settings devoid of prominent light or thermal sources, vibration energy harvesting emerges as a viable strategy [4][5][6][7][8][9].For the purpose of harvesting vibration energy, piezoelectric materials serve as an optimal medium for converting vibration energy into electrical power, efficiently harnessing minor mechanical energy that is frequently overlooked or wasted [10][11][12][13][14][15][16][17][18][19][20][21][22].Among virous piezoelectric transducers, the most commonly used structure in energy harvesting is the cantilevered beam or shaped cantilevered beam.A cantilevered piezoelectric energy harvester consists of a multilayer composite beam with one end secured to a host structure.In terms of layer arrangement, these harvesters primarily feature unimorph and bimorph configurations.
Numerous analytical models have been proposed to analyze the performance of cantilevered piezoelectric energy harvesters.Early efforts to model cantilevered piezoelectric energy harvesters drew inspiration from electromagnetic generator modeling [23,24], leading to the development of a lumped parameter model, which is often recognized as a single-degree-of-freedom model [25][26][27].While the lumped parameter model initially provided valuable insight into the electromechanical coupling performance of cantilevered piezoelectric energy harvesters, it oversimplified by reducing the system to a single vibration mode of the bender, neglecting key aspects like modal shape and strain distribution.Building upon previous work, Sodano et al. [28] introduced enhanced modeling techniques for piezoelectric energy harvesters, employing approximate distributed parameter models based on Hamilton's principle and the Rayleigh-Ritz method with Euler-Bernoulli beam assumptions.Following this path, Liao and Sodano further studied the optimization of cantilevered piezoelectric energy harvesters [29,30].Shortly afterward, Erturk and Inman identified limitations in the above models [31,32].Under Euler-Bernoulli beam assumptions, they then introduced fully coupled distributed parameter models with exact analytical solutions for these harvesters, incorporating both unimorph [33] and bimorph [34] configurations.Subsequently, considering Euler-Bernoulli, Rayleigh, and Timoshenko beam assumptions, Erturk developed a comprehensive framework for electromechanical modeling of piezoelectric energy harvesters utilizing the assumed-modes method [35].
In the aforementioned common modeling of cantilevered piezoelectric energy harvesters, the electric field strength is generally considered to be uniform, and it equals the output voltage divided by the piezoelectric layer thickness.This uniform electric field assumption seems reasonable at first glance, since the piezoelectric layer looks like a parallel-plate capacitor.However, for a piezoelectric bender, the strain distribution on a cross-section along the thickness direction is not uniform, which means the internal electric field generated by the spontaneous polarization should not be uniform.Following the electric field assumption of strongly coupled piezoelectric actuators [36], a non-uniform electric field is introduced to the modeling of strongly coupled piezoelectric bimorph harvesters by Gibus et al., resulting in a two-degree-of-freedom model [37].Later, driven by the electric field distribution derived using the exact solutions for plane problems [38], Wang and Shi pointed out the critical issues induced by the uniform electric field assumption and proposed a non-uniform electric field in the piezoelectric layer [39].In the research by Wang and Shi, an enhanced distributed parameter electromechanical model is developed, and the justification for incorporating non-uniform electric field consideration is systematically discussed.
With the introduction of a non-uniform electric field, improvements have been made to the conventional analysis models for piezoelectric energy harvesters.However, in current investigations, the performance correction arising from non-uniform electric field consideration is either neglected [37] or addressed only in specific instances involving certain sample harvesters with particular material parameters and geometric dimensions [39].This means that the general determinants of performance correction and the influence of accounting for non-uniform electric fields on the identification of material parameters remain unclear.
In the present study, adopting non-uniform electric field assumption in the piezoelectric layer, a unified distributed parameter electromechanical model considering both unimorph and bimorph configurations is presented and further simplified to a more practical single-mode one.Based on this, the general determinates of performance correction induced by non-uniform electric field consideration are proposed regardless of specific material parameters and geometric dimensions.As for an inverse problem, the influence of accounting for non-uniform electric fields on the identification of material parameters is also investigated.
The rest of the paper is outlined as follows: In Section 2, a unified distributed parameter model for cantilevered piezoelectric energy harvesters is introduced, which contains unimorph and bimorph configurations as well as uniform and non-uniform electric filed assumptions.In Section 3, the non-uniform electric field in the piezoelectric layer is validated by comparing the power output performance predicted using the present model with that obtained using two-dimensional finite element analysis.In Section 4, using the unified model presented in Section 2, the determinants of performance correction induced by incorporating non-uniform electric field assumption are revealed.In Section 5, as an inverse problem, the influence of accounting for non-uniform electric fields on the identification of material parameters is investigated.And Section 6 concludes the study.

Unified Modeling of Cantilevered Piezoelectric Energy Harvesters
Cantilevered piezoelectric energy harvesters are primarily composed of three configurations: the unimorph harvester (Figure 1a), the bimorph harvester with series connection (Figure 1b), and the bimorph harvester with parallel connection (Figure 1c).These configurations can be conceptualized as a composite beam, incorporating elastic and piezoelectric layers arranged in various configurations.
Sensors 2024, 24, x FOR PEER REVIEW practical single-mode one.Based on this, the general determinates of performance c tion induced by non-uniform electric field consideration are proposed regardless o cific material parameters and geometric dimensions.As for an inverse problem, the ence of accounting for non-uniform electric fields on the identification of material p eters is also investigated.
The rest of the paper is outlined as follows: In Section 2, a unified distributed p eter model for cantilevered piezoelectric energy harvesters is introduced, which co unimorph and bimorph configurations as well as uniform and non-uniform electri assumptions.In Section 3, the non-uniform electric field in the piezoelectric layer i dated by comparing the power output performance predicted using the present m with that obtained using two-dimensional finite element analysis.In Section 4, usin unified model presented in Section 2, the determinants of performance correction in by incorporating non-uniform electric field assumption are revealed.In Section 5, inverse problem, the influence of accounting for non-uniform electric fields on the i fication of material parameters is investigated.And Section 6 concludes the study.

Unified Modeling of Cantilevered Piezoelectric Energy Harvesters
Cantilevered piezoelectric energy harvesters are primarily composed of three c urations: the unimorph harvester (Figure 1a), the bimorph harvester with series co tion (Figure 1b), and the bimorph harvester with parallel connection (Figure 1c).configurations can be conceptualized as a composite beam, incorporating elastic an zoelectric layers arranged in various configurations.

Basic Equations
Considering the perspective of piezoelectricity, the constitutive equations for the elastic and the piezoelectric layers can be formulated as: Here, the subscript/superscript p and s indicate the piezoelectric and substructure layers, respectively.And the directions 1 and 3 align with the x and z directions.T 1 and S 1 denote stress and strain components along the x direction, respectively.Y represents the Young modulus, and Y p = 1/s E  11 , where s E 11 represents the compliance coefficient of the piezoelectric material under a constant electric field.d 31 represents the piezoelectric constant, and ε T 33 denotes the permittivity under constant stress.E 3 stands for the electric field, while D 3 stands for the electric displacement.In order to derive the one-dimensional material parameters s E 11 , d 31 , and ε T 33 from three-dimensional material parameters, it is necessary to employ a plane simplification approach and ensure zero stress along the 3-axis, as outlined in the Appendix A alongside the three-dimensional constitutive equations.
As an insulating material, the electric displacement must adhere to the electrostatic equilibrium equation, which can be expressed as: And the correlation between the potential distribution Φ and the electric field E 3 is expressed as: The strain component in the harvester could be expressed in terms of the radius of curvature according to Euler-Bernoulli beam assumptions as follows: Here, z denotes the position from the neutral axis, and w rel (x, t) represents the relative transverse deflection of the beam.When a harvester experiences a base motion, that reflected on the beam can be described as: where g(t) and h(t) are a traverse motion and a slight rotation at the base.Neglecting the extra excitation caused by the air, the mechanical governing equation of a cantilevered beam could be formulated as: In this expression, M(x, t) denotes the internal bending moment.w rel (x, t) and m represent the relative transverse deflection and the mass per unit length of the beam, respectively.c s I represents the equivalent damping term arising from structural viscoelasticity (where c s denotes the equivalent coefficient of strain rate damping and I is the equivalent area moment of the cross section), c a is the coefficient of viscous air damping, and these two distinct damping factors are assumed to account for the strain rate damping and viscous air damping, respectively.A comprehensive discussion on addressing the damping mechanism in cantilevered harvesters can be found in reference [32].The electrical governing equation of a composite piezoelectric beam subjected to base motion connected to an external resistance R l could be expressed as: where v(t) denotes the output voltage and q(t) represents the free charge generated by the energy harvester, which could be characterized by: Here, D represents the electric displacement vector on the surface of the piezoelectric layer, while n is the unit normal vector to the electrode.

Non-Uniform Electric Field in Cantilevered Piezoelectric Energy Harvesters
In traditional modeling of cantilevered piezoelectric energy harvesters, the electric field in the piezoelectric layer is generally considered to be uniform as the ratio of the output voltage to the thickness of the piezoelectric layer.However, if the electric displacement D 3 is calculated by substituting a uniform electric field E 3 and Equation (6) into Equation (2), the obtained electric displacement obviously does not adhere to the electrostatic equilibrium equation, i.e., Equation (4).Since the divergence of electric displacement signifies the density of free charges, the electrostatic equilibrium equation inherently implies that the free charge within an insulator must equate to zero.Consequently, within the Euler-Bernoulli beam framework, the assumption of a uniform electric field distribution contradicts a fundamental physical premise, namely, the absence of free charges within an insulator.
Herein, we combine Equations ( 1), ( 2), ( 4) and (6).A newly proposed non-uniform electric field can be resolved as: where Y p is the permittivity at constant strain The corresponding electric potential could be resolved as: where B(x, t) and C(x, t) are integration constants that can be solved using the potential boundary conditions at the surface of the piezoelectric layer, as illustrated in Figure 1.Specifically, for unimorph harvesters, the electric field and electric potential can be solved as: where h p is the piezoelectric layer thickness; h b and h c are the position of the bottom and top boundaries of the piezoelectric layer to the neutral axis, respectively; and h pc = (h b + h c )/2 is the position of the thickness center of the piezoelectric layer to the neutral axis.Defining the geometric ratio n g = h s /h p and elastic ratio n e = Y s /Y p , the positions of the layer boundaries of a unimorph harvester to the neutral axis, h a , h b , h c , and h pc , could be represented by n g , n e , and h p as: 2 n e +2n g n e +1 n g n e +1 2 n e +n g n e n g n e +1 Similarly, for bimorph harvesters with series connection, the electric field and electric potential could be solved as: For bimorph harvesters with parallel connection, the electric field and electric potential could be solved as:

Unified Electromechanical Governing Equations Containing Both Uniform and Non-Uniform Electric Field Assumption
Using the above basic equations combined with different assumptions of the electric field, i.e., uniform electric field and non-uniform electric field, the expression of the internal bending moment for both unimorph and bimorph configurations could be unified as [33,34,39]: This expression is formulated by calculating the internal bending moment of the beam using the integration of the normal stress on the cross-section, in which YI is the bending stiffness term for short-circuit condition and ϑ is the electromechanical coupling term.For a unimorph harvester, considering a uniform electric field, the bending stiffness term could be expressed using the width of the harvester b, the geometric and elastic parameters of the piezoelectric layer h p and Y p , the geometric ratio n g = h s /h p , and the elastic ratio n e = Y s /Y p : When employing the previously proposed non-uniform electric field, a refined bending stiffness can be derived as: As for the coupling term of a unimorph harvester, that obtained based on both uniform and non-uniform electric fields could be expressed as: For a bimorph harvester considering a uniform electric field, whether it is connected in series or parallel, the bending stiffness could be expressed as: When incorporating the non-uniform electric field proposed previously, the corrected bending stiffness of a bimorph harvester for both series and parallel connections could be derived as: As for the coupling term of a bimorph harvester with series connection, that obtained based on both uniform and non-uniform electric fields could be expressed as: And the coupling term of a parallel connected bimorph harvester is ϑ = −2ϑ B .It is noteworthy that if the grounded electrodes of a parallel connected bimorph harvester switch to the inner side, this coupling term would become 2ϑ B .For different electric field assumptions and different configurations, the expressions of YI and ϑ are summarized in Table 1.If the electrode and piezoelectric layer cover the entire length of the harvester, Equation ( 19) could be rewritten as where L is the length of the harvester, and H(x) is the Heaviside function.Substituting Equation ( 26) into Equation ( 8) yields the unified mechanical governing equation of motion with electric coupling: where δ(x) is the Dirac delta function, and it satisfies: For various electric field assumptions and different harvester configurations, the electrical governing equation could be unified by substituting the ungrounded side surface electric displacement, together with its outward unit normal vector, into Equation ( 9): where C p is the capacitance.For different configurations, C p is presented in Table 1.Based on the expressions for YI, ϑ, and C p shown in Table 1, it should be noticed that when the non-uniform electric field assumption is applied, only the bending stiffness term YI is affected, while the other parameters remain unchanged.

Dynamic Response to Harmonic Base Motion
To solve the electromechanically coupled equations given in Equations ( 27) and (29), the relative transverse deflection w rel (x, t) is represented by an absolutely and uniformly convergent series of the eigenfunctions, expressed as: where ϕ r (x) is the mass normalized eigenfunction, and η r (t) is the modal coordinate of the rth mode.For a clamped-free cantilevered beam with no tip mass, ϕ r (x) is provided as: where λ r is the dimensionless frequency numbers obtained from the characteristic equation given by: 1 and σ r is expressed as: Equation ( 32) can be numerically solved, and the first three orders of λ r are obtained as λ 1 = 1.8751, λ 2 = 4.6941, and λ 3 = 7.8548.
The mass normalized eigenfunction ϕ r (x) obeys the following orthogonality relations: Sensors 2024, 24, 4943 9 of 24 where δ rs is the Kronecker delta, equal to 1 when s = r and 0 otherwise, and ω r is the short-circuit undamped natural frequency of the rth mode, given by: Substituting Equation (30) into Equations ( 27) and ( 29), respectively, considering the orthogonality conditions mentioned in Equation ( 34), the electromechanical governing equations can be obtained as: where: χ r is the modal coupling term: and ζ r is the modal damping ratio, expressed as: The modal damping ratio ζ r is typically determined experimentally in practice, thereby eliminating the need to obtain the values of c s I and c a .A detailed discussion on resolving the damping coefficients from the modal damping ratio can be found in reference [33].
If the base motion w b (x, t) is harmonic in the form of w b (x, t) = (Z 0 + xθ 0 )e jωt (where j is the imaginary unit, ω denotes the driving circular frequency, and Z 0 and θ 0 are the amplitude of base translation and rotation, respectively), then the output voltage is also harmonic in the form of v(t) = V 0 e jωt (where V 0 is the amplitude of the output voltage).In this case, the electromechanical governing equations given by Equation (36) becomes: Solving Equation ( 40), the steady-state solution of the modal coordinate and the output voltage amplitude can be given as: And the relative transverse deflection w rel (x, t) can be easily obtained using Equation (30) and Equation (41).Although Z 0 and θ 0 are real-valued, v(t) and w b (x, t) do not have to be in phase, which makes V 0 a complex-valued quantity.Detailed development and resolutions of the electromechanical governing equations can be found in reference [39].

Reduced Single-Mode Response to Harmonic Base Translation
In most cases, the mode of interest is the fundamental vibration mode of the harvester (i.e., r = 1), which means that the excitation of the harvester is considered to be around ω 1 .Accordingly, by simplifying the base motion to a harmonic translation in the transverse direction (i.e., w b (x, t) = Z 0 e jωt ), the reduced expression for the amplitude of the voltage V * 0 across the external resistance can be written as: And the absolute value of V * 0 can be expressed as: Furthermore, the absolute value of power output can be expressed as: The formula of the output voltage and output power can be rewritten with respect to dimensionless variables as follows: where ω 2 Z 0 represents the amplitude of the base excitation acceleration, and ζ 1 , Ω, α, and k 2 are dimensionless variables.Specifically, ζ 1 is the damping ratio for the fundamental mode of the harvester, as presented in Equation (39).Ω is the frequency ratio, α is the resistance turning ratio, and k 2 is the squared effective electromechanical coupling coefficient.
Considering the mass normalized eigenfunction adopted in this study, Ω, α, and k 2 could be represented by the original system parameters as follows: Here, the effective electromechanical coupling coefficient squared k 2 is an important parameter directly reflecting the mechanical-to-electrical energy conversion efficiency of the system.k 2 describes the effectiveness of quasi-static energy conversion between electrical and mechanical forms.For a cantilevered energy harvester in an open circuit subjected to quasi-static stress, it is equal to the electrostatic energy divided by the total energy in the system.Theoretically, it is also related to the open-circuit undamped natural frequency (ω OC 1 ) and the short-circuit undamped natural frequency (ω SC 1 or ω 1 ) as: Sensors 2024, 24, 4943 11 of 24 Equation (49) provides a method to experimentally determine the squared effective electromechanical coupling coefficient k 2 while the damping ratio is far less than 1.0, approximately ς 1 < 0.2.
The optimal resistance turning ratio can be determined by differentiating the power expression Equation (47) with respect to the resistance tuning ratio, letting the derivative be zero, and solving the optimal resistance turning ratio: Substituting Equation (50) into Equation (47), the optimal power as a function of the frequency ratio can be obtained.The curve of the optimal power with respect to the frequency ratio or excitation frequency is considered as a power envelope.

Verification of the Non-Uniform Electric Field Assumption
In Section 2, considering both uniform and non-uniform electric field assumptions, the electromechanically coupled model of cantilevered piezoelectric energy harvesters is provided.While incorporating the non-uniform electric field, the bending stiffness term YI for both unimorph and bimorph harvesters is explicitly modified as presented in Table 1.To demonstrate that the model incorporating a non-uniform electric field in the piezoelectric layer is more accurate than the traditional model considering a uniform electric field, sample harvesters are employed and a numerical analysis is conducted using both the models, then compared with the results of a finite element analysis.
Herein, a unimorph sample harvester (H p = 0.5 mm, H s = 0.5 mm, b = 20 mm, L = 100 mm) and a parallel-connected bimorph sample harvester (H p = 0.5 mm, H s = 0, b = 20 mm, L = 100 mm) are introduced to demonstrate the performance corrections when incorporating a non-uniform electric field assumption.The sample harvesters are made from PZT-5H and aluminum.In Table 2, the one-dimensional material parameters of aluminum and PZT-5H under plane stress conditions are provided.Detailed resolutions of the one-dimensional material parameters of PZT-5H from three-dimensional material parameters for the plane problem are presented in the Appendix A. As for the damping ratio, we suppose ζ 1 = 0.01.In this section, a reduced single-mode response to harmonic base motion presented in Section 2.5 is employed, and the power output FRF (frequency response function) is defined as the output power P * 0 divided by squared excitation acceleration −ω 2 Z 0 2 .Power output prediction of these two sample harvesters based on a reduced single-mode model is conducted using MATLAB programs, which are provided in the Supplementary Materials.The finite element results for reference were achieved using COMCOL Multiphysics 5.6.In the COMSOL project, default PZT-5H material parameters, manually entered aluminum material parameters, a solid mechanics module, a 2D model, and plane stress analysis are employed.To adopt the Rayleigh damping mechanism, the first two orders of the damping ratio and the short-circuit undamped natural frequency are required.Therefore, we suppose ζ 1 = 0.01, ζ 2 = 0.02.And the first two orders of short-circuit undamped natural frequency vary with respect to specific sample harvesters, which can be obtained using steady-state analysis in the COMSOL project by setting the damping ratio and the circuit load resistance to zero.For the unimorph sample harvester in the COMSOL project, f 1 = 58.307Hz and f 2 = 365.24Hz.As for the bimorph sample harvester in the COMSOL project, f 1 = 46.952Hz and f 2 = 294.11Hz.
Numerical results of these two sample harvesters are presented in Figure 2. The first rows of Figure 2a,b represent the envelope of the power output FRF of these two sample harvesters with respect to the excitation frequency obtained based on Equations ( 47) and (50).The second row is the optimal resistance obtained based on Equation (50).And the last row is the power output FRF, while the load resistance is exactly the optimal resistance for the short-circuit undamped natural frequency ω 1 .Both the traditional model considering a uniform electric field as well as the improved model incorporating a non-uniform electric field are taken into account.The simulated power output results in COMSOL projects are provided in the last row for verification.

Determinants of Performance Corrections
According to the modeling of cantilevered piezoelectric energy harvesters incorporating non-uniform electric fields, the performance corrections to traditional model considering uniform electric fields are ultimately due to the bending stiffness correction, as presented in Table 1.As a key parameter, the bending stiffness term has a great influence on both the static and dynamic performance of cantilevered piezoelectric composites.For unimorph harvesters, the bending stiffness term   (Equation ( 20)) is modified to   ′ (Equation ( 21)).Similarly, for series-and parallel-connected bimorph harvesters, the bending stiffness term   (Equation ( 23)) is modified to   ′ (Equation ( 24)).The respective correction factors for the bending stiffness of unimorph and bimorph harvesters can be defined as: For the unimorph sample harvester, using the improved model incorporating a nonuniform electric field, the short-circuit undamped natural frequency is f 1 = ω 1 /(2π) = 58.289Hz, which is very close to the COMSOL result, and the corresponding optimal resistance is 11.28 kΩ.As for the parallel connected bimorph sample harvester, the shortcircuit undamped natural frequency is f 1 = 46.934Hz, which is also very close to the COMSOL result, and the corresponding optimal resistance is 4.10 kΩ.For the unimorph sample harvester connected to a 11.28 kΩ load resistance, the power output FRF is presented in the last row of Figure 2a, and that of the bimorph sample harvester connected to a 4.10 kΩ load resistance is presented in the last row of Figure 2b.Whether for a unimorph harvester or a bimorph harvester, it is obvious that the result obtained using the improved model incorporating a non-uniform electric field is basically consistent with the COMSOL 2D result, while a resonant frequency shift comparative to the traditional model is observed.

Determinants of Performance Corrections
According to the modeling of cantilevered piezoelectric energy harvesters incorporating non-uniform electric fields, the performance corrections to traditional model considering uniform electric fields are ultimately due to the bending stiffness correction, as presented in Table 1.As a key parameter, the bending stiffness term has a great influence on both the static and dynamic performance of cantilevered piezoelectric composites.For unimorph harvesters, the bending stiffness term YI U (Equation ( 20)) is modified to YI U ′ (Equation ( 21)).Similarly, for series-and parallel-connected bimorph harvesters, the bending stiffness term YI B (Equation ( 23)) is modified to YI B ′ (Equation ( 24)).The respective correction factors for the bending stiffness of unimorph and bimorph harvesters can be defined as: n g n e +1 n g 4 n e 2 +4n g 3 n e +6n g 2 n e +4n g n e +1 Recall that the electromechanical coupling coefficient of piezoelectric materials is defined as: Equation ( 51) can be rewritten as: n g n e +1 n g 4 n e 2 +4n g 3 n e +6n g 2 n e +4n g n e +1 Through Equation (53), it is found that the bending stiffness correction and the resulting performance correction are only related to three dimensionless parameters, i.e., the geometric ratio n g = h s /h p , the elastic ratio n e = Y s /Y p , and the electromechanical coupling coefficient squared k 31 2 .The impact of k  For a specific piezoelectric material, for example, PZT-5H, whose electromechanical coupling coefficient squared  31 2 = 0.151, the bending stiffness correction factor of unimorph harvester   with respect to geometric ratio   and elastic ratio   is presented For a specific piezoelectric material, for example, PZT-5H, whose electromechanical coupling coefficient squared k 31 2 = 0.151, the bending stiffness correction factor of unimorph harvester δ U with respect to geometric ratio n g and elastic ratio n e is presented in Figure 4.And the bending stiffness correction factor of bimorph harvester δ B with respect to n g and n e is presented in Figure 5.In both Figures 4 and 5, four discrete elastic ratios n e and continuous geometric ratios n g are considered, since the elastic parameter of the substrate is usually limited to a few commonly used materials, while the geometric dimensions have many variations.As we can see in Figures 4 and 5, it is found that both   and   are constantly greater than 1.0, and they decrease monotonically as   or   increases.While   is larger than 2.0, both   and   are less than 1.01.This means that, when the substructure  As we can see in Figures 4 and 5, it is found that both   and   are constantly greater than 1.0, and they decrease monotonically as   or   increases.While   is larger than 2.0, both   and   are less than 1.01.This means that, when the substructure As we can see in Figures 4 and 5, it is found that both δ U and δ B are constantly greater than 1.0, and they decrease monotonically as n g or n e increases.While n g is larger than 2.0, both δ U and δ B are less than 1.01.This means that, when the substructure layer's thickness is over twice the piezoelectric layer's thickness, the bending stiffness corrections for the PZT-5H unimorph and bimorph harvesters are less than 1%.This result is obtained when n e is in the range from 0.5 to 2.0, which covers most of the commonly used substrate materials.
On the other hand, for PZT-5H unimorph harvesters, while n g decreases to zero, a maximum bending stiffness correction factor δ U = 1.178 is observed, which indicates a 17.8% bending stiffness correction and an 8.5% short-circuit undamped natural frequency correction.For PZT-5H bimorph harvesters, while n g decreases to zero, a maximum bending stiffness correction factor δ B = 1.045 is observed, which indicates a 4.5% bending stiffness correction and a 2.2% short-circuit undamped natural frequency correction.From Figures 4 and 5, we know that a maximum bending stiffness correction factor can be achieved while the substructure layer thickness is zero.Then, substituting n g = 0 into Equation ( 53), the maximum bending stiffness correction factor is obtained as Through Equation (54), it is found that the achievable maximum bending stiffness correction and resulting maximum performance correction are only related to one dimensionless parameter, i.e., the electromechanical coupling coefficient squared, k layer's thickness is over twice the piezoelectric layer's thickness, the bending stiffness corrections for the PZT-5H unimorph and bimorph harvesters are less than 1%.This result is obtained when   is in the range from 0.5 to 2.0, which covers most of the commonly used substrate materials.On the other hand, for PZT-5H unimorph harvesters, while   decreases to zero, a maximum bending stiffness correction factor   = 1.178 is observed, which indicates a 17.8% bending stiffness correction and an 8.5% short-circuit undamped natural frequency correction.For PZT-5H bimorph harvesters, while   decreases to zero, a maximum bending stiffness correction factor   = 1.045 is observed, which indicates a 4.5% bending stiffness correction and a 2.2% short-circuit undamped natural frequency correction.From Figures 4 and 5, we know that a maximum bending stiffness correction factor can be achieved while the substructure layer thickness is zero.Then, substituting   = 0 into Equation (53), the maximum bending stiffness correction factor is obtained as ( ) Through Equation (54), it is found that the achievable maximum bending stiffness correction and resulting maximum performance correction are only related to one dimensionless parameter, i.e., the electromechanical coupling coefficient squared,  31 2 .The maximum bending stiffness correction factors    and    with respect to  31 2 are presented in Figure 6.Although the maximum correction factors for unimorph and bimorph harvesters are both achieved while the harvester consists of only piezoelectric material,    is obviously much larger than    at the same  31 2 .This is due to different internal electric field distributions.
0.00 0.02 0.04 0.06 0.08 0.10 0.12 0.14 0.16 0.18 0.20  Among commonly used piezoelectric materials, as pointed out in Figures 3 and 6, PZT-5H has the largest electromechanical coupling coefficient squared, i.e.,   Among commonly used piezoelectric materials, as pointed out in Figures 3 and 6, PZT-5H has the largest electromechanical coupling coefficient squared, i.e., k 063.This indicates a 25% maximum bending stiffness correction and a resulting 11.8% maximum short-circuit undamped natural frequency correction for unimorph harvesters, as well as a 6.3% maximum bending stiffness correction and a resulting 3.1% maximum short-circuit undamped natural frequency correction for bimorph harvesters.

Material Parameters Identification Considering Non-Uniform Electric Field
In Sections 3 and 4, the improved electromechanical model considering a non-uniform electric field in the piezoelectric layer is utilized to better predict the performance of cantilevered piezoelectric energy harvesters, in which case the material parameters are given in advance.Another scenario in which the improved model can make a difference is the identification of material parameters.
Case I: Only the Young modulus Y s of the substructure needs to be identified.In this case, the unimorph sample harvester introduced in Section 3 is employed (a PZT-5H and aluminum unimorph harvester: H p = 0.5mm, H s = 0.5mm, b = 20mm, L = 100mm).It is assumed that the material parameters of the piezoelectric layer, the geometric dimensions of the harvester, and the mass density of the substructure are known as listed in Table 2, and the Young modulus Y s of the substructure needs to be identified using the experimentally obtained overall characteristics of the energy harvester.The procedures to identify Y s are presented in Figure 7.In this case, the damping ratio ζ 1 and the short-circuit damped natural frequency f damped 1 are experimentally determined in advance.
Sensors 2024, 24, x FOR PEER REVIEW 18 of 26 maximum bending stiffness correction and a resulting 11.8% maximum short-circuit undamped natural frequency correction for unimorph harvesters, as well as a 6.3% maximum bending stiffness correction and a resulting 3.1% maximum short-circuit undamped natural frequency correction for bimorph harvesters.

Material Parameters Identification Considering Non-Uniform Electric Field
In Sections 3 and 4, the improved electromechanical model considering a non-uniform electric field in the piezoelectric layer is utilized to better predict the performance of cantilevered piezoelectric energy harvesters, in which case the material parameters are given in advance.Another scenario in which the improved model can make a difference is the identification of material parameters.
Case I: Only the Young modulus   of the substructure needs to be identified.In this case, the unimorph sample harvester introduced in Section 3 is employed (a PZT-5H and aluminum unimorph harvester:   = 0.5,   = 0.5,  = 20,  = 100).It is assumed that the material parameters of the piezoelectric layer, the geometric dimensions of the harvester, and the mass density of the substructure are known as listed in Table 2, and the Young modulus   of the substructure needs to be identified using the experimentally obtained overall characteristics of the energy harvester.The procedures to identify   are presented in Figure 7.In this case, the damping ratio  1 and the short-circuit damped natural frequency  1  are experimentally determined in advance.In order to demonstrate the accuracy of the material parameters identification considering a non-uniform electric field, a numerical experiment is conducted using the twodimensional COMSOL project, and the theoretically identified   is compared with the one originally set in the COMSOL project.
From the numerical experiment, the short-circuit undamped natural frequency of the unimorph sample harvester is obtained as  1 = 58.307Hz.For a real unimorph harvester, its damping ratio  1 and short-circuit damped natural frequency  1  can be obtained through a free vibration test, and the corresponding short-circuit undamped natural frequency could be obtained as If  1 is very small,  1 and  1 would be very close to each other.Combining the mass density and cross-section area of the unimorph harvester, its mass per unit length can be calculated as  = In order to demonstrate the accuracy of the material parameters identification considering a non-uniform electric field, a numerical experiment is conducted using the two-dimensional COMSOL project, and the theoretically identified Y s is compared with the one originally set in the COMSOL project.
From the numerical experiment, the short-circuit undamped natural frequency of the unimorph sample harvester is obtained as f 1 = 58.307Hz.For a real unimorph harvester, its damping ratio ζ 1 and short-circuit damped natural frequency f damped 1 can be obtained through a free vibration test, and the corresponding short-circuit undamped natural frequency could be obtained as  21), the elastic ratio can be solved as n e = 1.1565, and the Young modulus Y s of the substructure can be identified as Y s = n e Y p = 70.10GPa.Note that the original set Y s = 70 GPa, the identified Y s , and the corresponding error are presented in Table 3.According to Table 3, the error in material parameters identification was reduced from 4.51% to 0.14% when the non-uniform electric field was considered.Case II: The material parameters of the piezoelectric layer, Y p , d 31 , ε T 33 , need to be identified.This case is relatively complicated compared to Case I.In this case, the unimorph sample harvester introduced in the first case is also employed.It is assumed that the geometric dimensions, the material parameters of the substrate layer, and the mass density of the piezoelectric layer are known as listed in Table 2, and the material parameters of the piezoelectric layer Y p , d 31 , ε T 33 need to be identified using the experimentally obtained overall characteristics of the energy harvester.The procedures to identify Y p , d 31 , ε T 33 are presented in Figure 8.In this case, the damping ratio ζ 1 , the short-circuit damped natural frequency f damped 1 , and two sets of excitation frequencies f with corresponding power outputs FRF P * 0 / −ω 2 Z 0 2 need to be experimentally determined in advance.
Specifically, in the present study, these characteristics of the unimorph sample harvester are gathered through a numerical experiment conducted using the two-dimensional COMSOL project.And the following identified material parameters are compared with the ones originally set in the COMSOL project.Firstly, we need to identify the power output FRF curve of this unimorph sample harvester using its experimentally obtained overall characteristics.According to Equation (47), four parameters are required, i.e., f 1 , ζ 1 , k 2 , and α.In the numerical experiment, the short-circuit undamped natural frequency is determined as f 1 = 58.307Hz, and the damping ratio ζ 1 is considered as 0.01.These two parameters can be obtained through free vibration test for a real unimorph harvester, which was described in the first case.The resolution of the other two parameters, k 2 and α, needs two sets of excitation frequency and corresponding power output FRF to be substituted, respectively, into Equation (47).With the circuit load resistance R l = 11.28 kΩ, the power output FRF of this unimorph sample harvester with respect to the driving frequency is shown in the third row of Figure 2a as a set of blue dots.Take two sets of these blue dots, i.e., f = [58.5, 59.5] Hz, P * 0 / −ω 2 Z 0 2 = [2.090, 1.081] × 10 −4 Ws 4 /m 2 .The substitution of these two sets of data into Equation (47) generates two algebra equations, and the unknown parameters can be solved as k 2 = 0.0402 and α = 0.4484.Recall the damping ratio ζ 1 = 0.01 and the short-circuit undamped natural frequency ω 1 = 58.307× 2π rad/s.The power output FRF of this unimorph sample harvester with respect to the excitation frequency can be identified using Equation (47) and is presented in Figure 9.
damped natural frequency  1  , and two sets of excitation frequencies  with corresponding power outputs FRF | 0 * | (− 2  0 ) 2 ⁄ need to be experimentally determined in advance.Specifically, in the present study, these characteristics of the unimorph sample harvester are gathered through a numerical experiment conducted using the two-dimensional COMSOL project.And the following identified material parameters are compared with the ones originally set in the COMSOL project.Firstly, we need to identify the power output FRF curve of this unimorph sample harvester using its experimentally obtained overall characteristics.According to Equation (47), four parameters are required, i.e.,  1 ,  1 ,  2 , and .In the numerical experiment, the short-circuit undamped natural frequency is determined as  1 = 58.307Hz , and the damping ratio  1 is considered as 0.01.These two parameters can be obtained through free vibration test for a real unimorph harvester, which was described in the first case.The resolution of the other two parameters,  2 and , needs two sets of excitation frequency and corresponding power output FRF to be substituted, respectively, into Equation (47).With the circuit load resistance   = 11.28 kΩ, the power output FRF of this unimorph sample harvester with respect to the driving frequency is shown in the third row of Figure 2a as a set of blue dots.Take two sets of these blue dots, i.e.,  = [58.⁄ .The substitution of these two sets of data into Equation (47) generates two algebra equations, and the unknown parameters can be solved as  2 = 0.0402 and  = 0.4484 .Recall the damping ratio  1 = 0.01 and the short-circuit undamped natural frequency  1 = 58.307× 2 rad s ⁄ .The power output FRF of this unimorph sample harvester with respect to the excitation frequency can be identified using Equation (47) and is presented in Figure 9. Since the power output FRF curve has been successfully identified as shown in Figure 9, the determined value of the resistance turning ratio  and the effective electromechanical coupling coefficient squared  2 should be trustworthy, i.e.,  = 0.4484 and  2 = 0.0402 .Furthermore, the circuit load resistance   = 11.28 kΩ , the short-circuit undamped natural circular frequency  1 = 58.307× 2 rad s ⁄ , the capacitance   , and the coupling term  of this unimorph harvester can be solved using Equation (48), which gives   = 1.0850 × 10 −7 F and  = 8.8757 × 10 −5 Nm V ⁄ .Recall that the bending stiffness was identified in the first case as  = 0.1107 N • m 2 .The expressions of , , and   actually formed three algebra equations containing three unknown constants, i.e.,   ,  31 , and  33  .While the traditional uniform electric field is adopted, these three algebra equations can be summarized as: Since the power output FRF curve has been successfully identified as shown in Figure 9, the determined value of the resistance turning ratio α and the effective electromechanical coupling coefficient squared k 2 should be trustworthy, i.e., α = 0.4484 and k 2 = 0.0402.Furthermore, the circuit load resistance R l = 11.28 kΩ, the short-circuit undamped natural circular frequency ω 1 = 58.307× 2π rad/s, the capacitance C p , and the coupling term ϑ of this unimorph harvester can be solved using Equation (48), which gives C p = 1.0850 × 10 −7 F and ϑ = 8.8757 × 10 −5 Nm/V.Recall that the bending stiffness was identified in the first case as YI = 0.1107 N•m 2 .The expressions of YI, ϑ, and C p actually formed three algebra equations containing three unknown constants, i.e., n e , d 31 , and ε T 33 .While the traditional uniform electric field is adopted, these three algebra equations can be summarized as: The unknown constants can be solved as: And the Young modulus of the piezoelectric layer can be obtained as Y p = Y s /n e = 63.11GPa.While the present non-uniform electric field is incorporated, these three algebra equations can be summarized as: The unknown constants can be solved as: In summary, the material parameters originally set in the numerical experiment and the identified material parameters are presented and compared in Table 4.As we can see in Table 4, while the non-uniform electric field is incorporated, the identification error of the permittivity at constant stress ε T 33 is hardly affected.However, the identification error of the Young modulus Y p is reduced from 4.12% to 0.30%, and that of the piezoelectric constant d 31 is reduced from 2.39% to 0.44%.The above procedures to identify material parameters are realized using MATLAB programs, which are provided in the Supplementary Materials.

Conclusions
Considering both unimorph and bimorph harvesters, an electromechanically coupled model incorporating a non-uniform electric field in the piezoelectric layer is presented and simplified to a more practical single-mode model.The present model is validated by comparing its power output prediction with that obtained using the COMSOL project.The results show that, when the non-uniform electric field in the piezoelectric layer is incorporated, the accuracy of both the performance prediction and the material parameters identification for cantilevered piezoelectric energy harvesters is improved.
The above constitutive equations can be simplified to two-dimensional ones for the plane problem.Specifically for the plane strain problem, substituting S 2 = S 23 = S 12 = 0 and E 2 = 0 into Equations (A1) and (A2) yields:  For the Euler-Bernoulli beam model adopted in the present study, an additional assumption T 3 = 0 can be applied and substituted into Equations (A3) and (A4), which lead to the one-dimensional constitutive equations of piezoelectric materials for the plane strain and plane stress problems, respectively: Sensors 2024, 24, 4943 23 of 24 For the typical PZT-5H material employed in the sample harvesters: In the main text, the one-dimensional material parameters for the plane stress problem, s E 11 , d 31 , and ε T 33 , are employed, and the values of these parameters are listed in Table 2.
Bimorph harvester with parallel connection.

Figure 2 .
Figure 2. Comparison of power output FRF for sample harvesters.

Figure 2 .
Figure 2. Comparison of power output FRF for sample harvesters.

Figure 3 .
Figure 3. Impact of material electromechanical coupling strength.

Figure 3 .
Figure 3. Impact of material electromechanical coupling strength.

Figure 4 .
Figure 4. Bending stiffness correction factor δ U with respect to geometric ratio n g and elastic ratio n e (PZT-5H unimorph harvesters).

Figure 5 .
Figure 5. Bending stiffness correction factor δ B with respect to geometric ratio n g and elastic ratio n e (PZT-5H bimorph harvesters).

6 .at the same k 31 2 .
Although the maximum correction factors for unimorph and bimorph harvesters are both achieved while the harvester consists of only piezoelectric material, δ Max U is obviously much larger than δ Max B This is due to different internal electric field distributions.Sensors 2024, 24, x FOR PEER REVIEW 17 of 26

Figure 6 .
Figure 6.Maximum bending stiffness correction factors    and    with respect to  31 2 .
31 2 = 0.151, and the corresponding maximum bending stiffness correction factors are as follows:    = 1.178 and    = 1.045.If  31 2 increases to 0.2, which is very likely for PZTbased composites such as macro-fiber composites, the maximum bending stiffness correction factors would become    = 1.250 and    = 1.063 .This indicates a 25%

Figure 6 .
Figure 6.Maximum bending stiffness correction factors δ Max U and δ Max B with respect to k 31 2 .

31 2 == 1 .
0.151, and the corresponding maximum bending stiffness correction factors are as follows: δ Max U = 1.178 and δ Max B = 1.045.If k 31 2 increases to 0.2, which is very likely for PZT-based composites such as macro-fiber composites, the maximum bending stiffness correction factors would become δ Max U 250 and δ Max B = 1.

Figure 7 .
Figure 7. Procedures to identify the Young modulus   of the substructure.

Figure 7 .
Figure 7. Procedures to identify the Young modulus Y s of the substructure.

Figure 9 .
Figure 9. Numerical experiment result and identified curve.

Figure 9 .
Figure 9. Numerical experiment result and identified curve.

Table 1 .
Expressions of YI, ϑ, and C p considering multiple configurations and different electric field assumptions.

Table 2 .
One-dimensional material parameters of PZT-5H and aluminum under plane stress conditions.
Through Equation (53), it is found that the bending stiffness correction and the resulting performance correction are only related to three dimensionless parameters, i.e., the geometric ratio   = ℎ  ℎ  ⁄ , the elastic ratio   =     ⁄ , and the electromechanical coupling coefficient squared  312.The impact of  31 2 on the bending stiffness correction factor   and   are reflected by the same term  31 31  2on the bending stiffness correction factor δ U and δ B are reflected by the same term k Sensors 2024, 24, x FOR PEERREVIEW  16 of 26in Figure4.And the bending stiffness correction factor of bimorph harvester   with respect to   and   is presented in Figure5.In both Figures4 and 5, four discrete elastic ratios   and continuous geometric ratios   are considered, since the elastic parameter of the substrate is usually limited to a few commonly used materials, while the geometric dimensions have many variations.
Sensors 2024, 24, x FOR PEERREVIEW  16 of 26in Figure4.And the bending stiffness correction factor of bimorph harvester   with respect to   and   is presented in Figure5.In both Figures4 and 5, four discrete elastic ratios   and continuous geometric ratios   are considered, since the elastic parameter of the substrate is usually limited to a few commonly used materials, while the geometric dimensions have many variations.
(20)is very small, f 1 and f damped 1 would be very close to each other.Combining the mass density and cross-section area of the unimorph harvester, its mass per unit length can be calculated as m = 0.102 kg/m.Substituting ω 1 = 58.307×2πrad/s,L=0.1 m and m = 0.102 kg/m into Equation(35), the bending stiffness YI can be solved as 0.1107 N•m 2 .Based on a traditional uniform electric field, substituting Y p , b, h p , n g , YI into Equation(20), the elastic ratio can be solved as n e = 1.2071, and the Young modulus Y s of the substructure can be identified as Y s = n e Y p = 73.16GPa.Based on a present non-uniform electric field, substituting Y p , b, h p , n g , d 31 , ε S 33 , YI into Equation (

Table 3 .
Identified Y s and corresponding error.

Experimentally obtained overall characteristics Two set of excitation frequency and corresponding power output FRF Short-circuit undamped natural frequency Solve: Resistance turning ratio and squared effective electromechanical coupling coefficient Using Eq.(47) Solve: Bending stiffness Using Eq.(35) Solve: Coupling term and capacitance Short-circuit damped natural frequency Damping ratio Using Eq.(48) Identify power output FRF curve and compare with experimental results Using Eq.(47) Using Eq.(20) (22) and expression of Using Eq.(21) (22) and expression of Solve: , and Uniform electric field: Non-uniform electric field: Solve: , and Figure 8.
Procedures to identify   ,  31 , and  33  .Figure 8. Procedures to identify Y p , d 31 , and ε T 33 .
bh p