Piezoelectric elements subjected to low frequency excitation. Empirical determination of stress and frequency influence on piezoelectric parameters

Designing optimum energy harvesting devices is the aim of several developments based on numerical or analytical studies of different piezoelectric configurations that usually consider constant piezoelectric properties. Experimental tests on bending piezoelectric patches showed that the electrical response depended on the frequency and amplitude of the mechanical excitation for displacement-imposed systems. Analytical and numerical calculations required adapting piezoelectric parameters to properly represent experimental results. A novel formulation to calculate piezoelectric parameters using the mechanical stress and the excitation frequency as inputs is proposed and discussed. A linear dependency on the mechanical stress of the piezoelectric ceramic and a logarithmic dependency on the excitation frequency have been combined to propose a unique calculation procedure. Later, this procedure was applied to compute different piezoelectric parameters to set numerical (2% error) and analytical (1% error) calculations that accurately represented experimental results. Finally, the practical implications of considering or not considering the frequency and stress dependency of the piezoelectric properties was evaluated for a theoretical bimorph cantilever configuration, whose excitation frequency decreased whereas the amplitude was kept constant. Results showed that only 1/3 of the energy production that was predicted with constant piezoelectric properties can be expected when considering frequency and stress influence.


Introduction
Piezoelectric energy harvesting is one of the most prolific research fields in the recent years. It has moved from 94 publications in 2010 to 948 publications in 2019 only considering Elsevier publication lists containing both concepts ''energy harvesting '' and ''piezoelectric''. Direct piezoelectric effect has been deeply explored in sensor development field. This application has turn into health and bio-chemical sensors (Platt et al. 2005;Park et al. 2008;Lu et al. 2013;Wu et al. 2019) in the recent years. Complementary, the faster development of nanotechnology applications, whose power supply requirements are in the range of mW or lW, has allowed using piezoelectric elements as energy suppliers. Examples of these low-consumption autonomous devices are micropumps (Ma et al. 2008), internal drug delivery systems (Staples et al. 2006), self-powered strain sensors (Huo et al. 2020), or asphalt self-powered temperature sensors (Hwang et al. 2019).
Two approaches gather most of the recent developments on piezoelectric energy harvesting: piezoelectric ceramic disks subjected to mechanical compressive efforts and piezoelectric ceramic patches subjected to bending efforts. The compressive approach is typically produced by piling piezoelectric disks , resulting in specific design options depending on the particular application, e.g. pavement energy harvesting . The compressive configuration is characterised by its elevated stiffness (Zhao et al. 2012) that promotes stress concentration problems that can be overcome by centring load devices (Guo and Lu 2019). Regarding the bending configuration, most of the proposed systems (Mitcheson et al. 2004; Moon et al. 2014;Yang et al. 2017a) are based on placing a vibrational mass at the end of a cantilever beam with piezoelectric patches installed. Combining bending and compressive piezoelectric configurations (Pérez-Lepe et al. 2016) led to Mooney, Cymbal or Bridge configurations. Among them, some authors pointed out that bridge configuration is the most effective one (Zhao et al. 2012). Possible improvements of these last energy harvesting systems are piling bridge devices (Jasim et al. 2018) or designing asymmetric cymbals (Goh et al. 2017).
Piezoelectric properties of a wide range of piezoelectric materials have been characterised. Caliò et al. (2014) comprehensively summarised part this information in a significative work. Many other researchers focused on studying specific properties, like the electro-mechanical coupling coefficient (Shu and Lien 2006), or the properties of specific piezoelectric materials like PMN-PT (Pramanik et al. 2019) or composites (Banerjee et al. 2015). However, only a few researchers have pointed out the idea that these properties are not constant and they depended on the frequency of the mechanical excitation (Damjanovic 1997;Fernandes et al. 2002) or the mechanical stress (Gusarov et al. 2016). Thus, an adaptative calibration of the mechanical-electrical response may be required for realistic design of piezoelectric energy harvesting devices.
The optimisation of piezoelectric energy harvesting devices tend to search for mechanical designs that amplify the excitation frequency (Umeda et al. 1996) or the mechanical stress on the piezoelectric elements (Saxena et al. 2017;Yildirim et al. 2017). Although these are intrinsically correct approaches, not including the influence of these parameters in the calculation of piezoelectric properties may cause deviations on the assessment of the effectiveness of the system respect to the theoretical designs.
In general, numerical simulations of piezoelectric devices did not consider the influence of the excitation frequency or the mechanical stress level on piezoelectric parameters, as per authors knowledge. In this line, Guo et al. (2019) andJasim et al. (2018) compared laboratory tests with simulations of road energy harvester systems. Toyabur et al. (2017) simulated cantilever systems with multiple degrees of freedom and Nowak et al. (2020) numerically studied the influence of design parameters of a bimorph cantilever system.
Most of the existing analytical models do not consider the influence of the mechanical stress and the frequency of the mechanical excitation on the definition of the piezoelectric parameters either. Eggborn (2003) combined Euler-Bernoulli's beam theory with fundamental piezoelectric definitions to obtain a simple analytical model for predicting voltage output. Yang et al. (2017b) focused on the analytical modelling of the energy conversion and Townley (2009) aimed to properly model experimental vibrational tests on cantilever generators. In the recent years researchers focused on the analytical modelling of the connection of piezoelectric elements (Basutkar 2019) or the study of composite configurations (Keshmiri et al. 2019).
According with the literature review there are a lot of studies proposing different energy harvesting designs (Caliò et al. 2014) but only a few studies pointing out that piezoelectric properties depend on the frequency of the mechanical excitation and mechanical stress in piezoelectric materials (Damjanovic 1997;Fernandes et al. 2002;Gusarov et al. 2016). There is no previous publication, as per authors knowledge, that combined the dependency of the piezoelectric properties on the excitation frequency and the mechanical stress in a single formulation as it is intended in this work. Hence, the research presented herein was aimed: (i) to provide additional experimental evidences on this dependency; (ii) to propose a novel and easy to implement formulation to represent it and (iii) to use the proposed formulation to foresee the influence of considering or not considering the influence of stress and excitation frequency on the evaluation of a theoretical design of an energy harvesting device. Experimental tests, numerical simulations and analytical calculations were carried out.

Materials and specimens
Experimental tests were conducted on ''cantilever'' unimorph specimens cyclically excited at a fixed frequency (three different frequencies per test configuration) and amplitude (20 mm). This amplitude was translated into different stress levels because of considering different lengths (test configuration) of the beams where the piezoelectric patches were mounted on. Thus, cantilever free vibration was not considered but motion was externally imposed instead to keep the peak stress values constant during tests. Two substrate pieces were used: an aluminium 6082 plate with a rectangular transversal section of 40 mm width and 2 mm thickness; and a high impact polystyrene plate with the same sectional dimensions. Young's modulus of polystyrene was experimentally determined by tensile tests (three repetitions) resulting to be 1.5 GPa. Young's modulus of aluminium plates was 68.9 GPa according with producer's datasheet.
Cantilever configuration was considered because of its clear boundary conditions and because it is easy to be modelled with analytical tools. Those facts allowed to focus the research on the study of the influence of the frequency and stress of the mechanical excitation on the piezoelectric parameters.
Piezoelectric P-876.A12 patches were purchased. This patches used PIC255 piezoelectric ceramic, which was made of modified lead zirconate-lead titanate and was classified type 200 (Soft PZT) according with EN 50324-1 (European Committee for Electrotechnical Standarization 2002). Patches and ceramic details are included in Table 1. Extended information about properties or production procedures may be accessed in (Physik Instrumente (PI) GmbH & Co. 2008. A piezoelectric patch P-876.A12 was bonded with cyanoacrylate to the substrate plate at 50 mm from one of the endings to leave the required free space in order to properly restrain the movement of this end of the plate by clamping it. Pressure and constant temperature were maintained during the curing of this adhesive (4 h). The restrained end of the plate was clamped to a fixed support. The connection state was checked before and after every test to assure cantilever configuration. Two length dimensions (free length was measured between the restrained edge and the displacement application point) were considered for aluminium plate (250 mm and 350 mm, corresponding to free end rotation angles of 0.120 and 0.086 respectively) and three for polystyrene plate (250 mm, 350 mm and 450 mm, corresponding to free end rotation angles of 0.120, 0.086 and 0.067 respectively). Tests were carried out at indoor environment with constant temperature (21°C ± 2°C) and relative humidity (65% ± 5%).

Experimental tests
The displacement at the ''free'' end of the plates was imposed by the eccentric rotation movement provided by an electric motor, whose rotation speed could be regulated. This rotation was transformed into a vertical linear movement (20 mm amplitude for all tests) by a steel tool (Fig. 1, left). Horizontal movement of the ''free'' end was restrained by vertical aluminium profiles at both sides of the plate (Fig. 1 right). Thus, the motion at the extreme of the beam was imposed to be a sinusoidal displacement of 20 mm amplitude. The corresponding frequency was regulated. The mechanical device of the testing setup is shown in Fig. 1, top and the testing configuration scheme is shown in Fig. 1, bottom. All tests were repeated for 4 different electrical loads (120 kX, 230 kX, 402 kX and 570 kX) to empirically select the optimum one (maximum energy generation). The optimum electrical load was set to 230 kX for the combination of tested frequencies (below 10 Hz) and amplitude (20 mm). This load was over the typical value assumed to represent the open circuit condition: 1/wC = 177 kX, calculated for the maximum excitation frequency initially planned (w = 10 Hz = 20p rad/s) and the capacitance of the piezoelectric (C = 90nF) reported by the manufacturer. The results of the tests with electrical loads that were different from 230 kX are not discussed in this document. All of them showed the same qualitatively response than the one described in detail for 230 kX load, but generated less energy output, so lower resolution of the analysed dependency.
This optimum loading resistance was set by installing two electric resistances, 220 kX and 10 kX, in series. This configuration (voltage divider) was set to limit the input voltage range into the data acquisition system at ± 10 V by measuring the electric signal on the smaller resistance (10 kX). Piezoelectric voltage output was continuously recorded at 200 Hz using a general-purpose data acquisition system (HBM Spider 8).
Mechanical excitation frequency was limited to 10 Hz because tests were part of a larger research project aimed to develop a new energy harvesting device to be installed in roads. The frequency of vehicles passing on it was calculated to be in the range of 1-2 Hz. In addition, this particular application is characterised by the fact that the resonance frequency of the system is far greater than the mechanical excitation frequency. Nevertheless, cantilever beam did not vibrate freely but an imposed and controlled oscillation was externally applied. Precise mechanical frequency excitation was obtained from voltage waveforms output signal for every testing case. Table 2 summarises the tests carried out. Test name, plate material, free cantilever length, excitation frequency and output voltage amplitude are included. Test name is in the form X_Y_Z, where X represents the plate material (A for aluminium and P for polystyrene), Y represents the cantilever free length (250 mm, 350 mm or 450 mm) and Z represents the different testing frequency (F1-F3 from the lower value to the greater one).

Numerical simulations
Numerical simulations were carried out using commercial software ANSYS 19.2 with MEMS Add-in. The geometric definition included the substrate plate and the piezoelectric patch, distinguishing the foil cover (Dupont Kapton Polyimide Film) from the ceramic material in the patch definition. That allowed to consider the realistic volume of the piezoelectric ceramic placed in the correct position. Dimensions of the piezoelectric patch and the piezoelectric ceramic were provided by the manufacturer. A three-dimensional model was used.
Linear elastic response was considered for the mechanical simulation of all materials. Ceramic volume was modelled with a simplified piezoelectric body whose properties were provided by the manufacturer and also included in some researches (Krommer et al. 2012; Physik Instrumente (PI) GmbH & Co. 2016). All material properties used in numerical simulations are summarised in Table 3. e31, e33, e15, ep11 and ep33 were only used in numerical simulations. g31 was the only parameter used in the proposed simplified analytical calculations.
Contacts between every pair of parts (substrate-foil and foil-piezoelectric ceramic) were assumed to be completely bonded, so no separation, sliding neither penetration were allowed. Mesh was composed of hexahedral (SOLID186) and tetrahedral (SOLID187) elements to discretize the substrate (5 mm size), hexahedral elements (SOLID226) for the piezoelectric ceramic (0.5 mm size) and tetrahedral elements (SOLID187) for the foil part around the ceramic (1 mm size). Convergence analysis on the mesh size was performed to set the sizing. A variation of strain below 5% between the used size and the following refinement (half size) was checked. Meshes for 250 mm, 350 mm and 450 mm plate cases had 30,670, 30,830 and 30,990 elements respectively.
The mechanical boundary condition was set by fixing the top and bottom faces of the plate substrate in a length of 50 mm to accurately represent experimental tests. Load was applied on the opposite edge as an imposed displacement orthogonal to the substrate plate plane with a value of 20 mm. Voltage of the bottom face of the piezoelectric ceramic (the closest to the substrate) face was set to 0 V and the voltage of the top face of the piezoelectric ceramic (the furthest from the substrate face) was coupled in all nodes. The voltage value on this face is the output result considered for discussion.
Static structural analysis was performed on the defined model. Calculation time was between 1 0 and 2 0 using an Intel Ò Core TM i7-9700 K CPU @3.60 GHz and 16.0 GB RAM memory (Fig. 2).

Analytical calculations
Analytical calculations were performed on the hypothesis of: (i) Euler-Bernoulli's beam theory, (ii) the well-known concepts of hybrid beams subjected to bending efforts and (iii) the linear relationship between the average axial stress in the piezoelectric ceramic and the output voltage that is defined through the piezoelectric voltage coefficient (g31). The following procedure was implemented: 1. Defining the ratio of the Young's modulus of the foil and Young's modulus of the piezoelectric ceramic respect to the substrate one (n f , n p respectively). 2. Calculating the equivalent width of foil and the equivalent width of the piezoelectric ceramic by homogenising the section to the substrate material.
Data provided by piezoelectric manufacturer 3. Calculating mechanical properties (area and moment of inertia) of the equivalent homogenised section. 4. Calculating the force that was necessary to cause the imposed displacement considering the cantilever configuration and the equivalent mechanical properties. The corresponding bending moment was also calculated. 5. Calculating the mechanical axial stress in the piezoelectric assuming a linear stress distribution on the homogenised section. 6. Calculating the output voltage considering the mechanical stress and the piezoelectric voltage coefficient (g31).

Results and discussion
Previously described procedures and properties are used to obtain experimental, numerical and analytical results. These are presented and discussed in this section. Next, piezoelectric parameters are adjusted to fit the numerical and the analytical calculations to the experimental results. This empirical fitting is analysed on the basis of the existing literature about the dependence of piezoelectric properties on mechanical excitation frequency and mechanical stress. Finally, a theoretical design of a piezoelectric energy harvesting device is evaluated with the initial piezoelectric parameters indicated by the provider and with the ones calculated considering frequency and stress dependency. This comparison allowed to quantify the influence of considering these effects of piezoelectric properties on energy harvesting.

Results
The frequency of the mechanical excitation and the output voltage amplitude are provided in the last two columns of Table 2. Table 4 summarises the numerical (FEA) and analytical results and it provides the corresponding relative errors respect to experimental results. These first calculations did not consider the influence of the frequency of the mechanical excitation or the mechanical stress. Experimental results are plotted in Fig. 3. It was observed that the piezoelectric output voltage depended on the excitation frequency for all tests. Voltage output increased when excitation frequency did, even for a fixed displacement amplitude and maintaining the same vibration shape, which was forced by the imposed external movement. Thus, it was expected that the calculation approaches that did not consider the influence of the frequency of the mechanical excitation of the piezoelectric patch did not represent their electrical response properly. It was the case of both numerical and analytical methods previously proposed, whose voltage output prediction was unique per each test configuration. Relative error moved from 22 to 81% in the case of numerical simulations, which tended to overestimate the output voltage. Analytical calculations tended to underestimate the output voltage in a range between 18 and 392%.
Both numerical and analytical calculation procedures showed their voltage output results were influenced by the substrate plate length in the same sense than experimental results were but with different magnitude. This fact proved that mechanical stress calculation of the piezoelectric ceramic, which was independent from the excitation frequency, was coherent but mechanical-electrical coupling piezoelectric properties were not.

Fitting calculation parameters
A four-steps procedure was followed to adjust the piezoelectric input parameters. This input parameters were different variables for analytical than for finite element calculations. First, the best-fitting value of these parameters was obtained for every study case. Second, the tendency of the value of the considered parameter was fitted with a suitable formulation for each substrate type (aluminium or polystyrene). Third, this formulation was used to calculate the corresponding input parameters for all studied cases. Finally, calculations of the output voltage were repeated considering these new set of parameters and the new results were analysed. This procedure was implemented twice: for numerical and for analytical analyse because different input parameters were required.

Fitting input parameters for numerical simulations
A previous initial step was added to the procedure described in Sect. 3.2 for the numerical simulations approach. The influence of the different input piezoelectric parameters (e31, e33, e15, ep11 and ep33) was analysed through a numerical sensitivity analysis. Results showed that modifying the values of e33, e15 and ep11 (100%, 300% and 40% respectively) had little influence (5%, 0% and 0% respectively) on the simulated output voltage. e31 and ep33 showed significant influence for the used piezoelectric configuration and both them were adjusted.
An additional numerical simulation was implemented considering e31 = -3C/m 2 and ep33 = 2042, which underestimated the output voltage and it was used as the second point for the linear interpolation. A pair of values for the parameters e i 31 and ep i 33 were obtained by linear interpolation between the initial numerical results (overestimated output voltage) and these additional numerical simulations (underestimated output voltage), to fit to the experimental results for every test. Interpolated optimum e i 31 and ep i 33 are summarised in 4th and 5th columns of Table 5.
According with (Damjanovic 1997), d33 piezoelectric coefficient showed a linear relationship with the logarithm of the mechanical excitation frequency. In addition, (Gusarov et al. 2016) proved that g31 piezoelectric voltage coefficient linearly depended on the mechanical stress of the piezoelectric ceramic. The following formulation is newly proposed for the calculation of e31 piezoelectric coefficient by combining both previous literature approaches: The ratio between e31 value provided by the piezoelectric patch manufacturer and the adjusted one, e a 31, was used to calculate the adjusted value of ep a 33 imposing the inverse proportion: The plot of e i 31 vs. ln(x) (Fig. 4 is an example of the procedure) was used to fit a linear equation for each specific plate (particular length and substrate material combination), resulting in k 4 e * (slope) and k 3 e * (independent term) values summarised in Table 6. Then the linear equations for the same substrate plate material (but different lengths) were modified by imposing a unique independent term per material, k 3 e . This independent term was obtained as the average of the independent terms of the previous linear fitting equations for the same substrate (aluminium or polystyrene). Slope coefficient, k 4 , was modified accordingly. At this point, different plate lengths of the same material had the same k 3 e coefficient but different k 4 e coefficients. Then, k 4 e coefficients were plot against r r N dimensionless value (Fig. 5 is an example of the procedure) where r was the average mechanical stress in the piezoelectric ceramic in the poling direction (this is a result of the first numerical simulations) and r N-= 44.2 MPa is the uniformly distributed equivalent stress associated to the blocking force (265 N for P-876.A12 piezoelectric patch). Linear fitting of k 4 e vs. r r N allowed to obtain k 1 e and k 2 e values completing the definition of the terms in Eq. 1. With the k 1 e , k 2 e and k 3 e parameters adjusted (see Table 6), piezoelectric coefficients, e a 31 and ep a 33 were calculated using Eqs. 1 and 2 for every specific testing case (see 6th and 7th columns in Table 5). Simulations were run again considering these piezoelectric parameters, which were calculated as a function of the stress in the piezoelectric ceramic (r) and the frequency of mechanical excitation. Output voltages are summarised in the 8th column of Table 5, V SIM . The average error of these simulations in terms of output voltage was -2.3% for aluminium plates and -1.8% for polystyrene plates (last column in Table 5). The average error (with no sign) was 6.5%. It was significantly lower than the value of 143.3% that was obtained with the numerical simulations that considered producer's properties. Thus, it was proved that the proposed methodology for adjusting piezoelectric properties as function of mechanical excitation frequency and stress level increased the accuracy of the numerical simulations.

Fitting input parameter for analytical calculations
Results of the analytical calculations performed using manufacturer's piezoelectric g31 parameter are summarised in the 4th column of Table 7. These values showed an average relative error (with no sign) of 91.3%. The optimum value of g31 parameter, g i 31, was interpolated for every case on the basis that the output voltage is proportional to it. Then, the same procedure that was described to compute k 4 e *, k 3 e *, k 4 e , k 3 e , k 2 e and k 1 e coefficients for e a 31 calculation was followed but to calculate g a 31 parameter according with Eq. 3. The same dependency on mechanical excitation frequency and mechanical stress level was set. Results of the fitting procedure (k 1 g , k 2 g and k 3 g coefficients to calculate g a 31 for every tested case) and the output voltage of the adjusted numerical model, V Ana , are summarised Tables 7 and 8.
The average relative error (with no sign) was reduced from 91.3% (obtained using manufacturer's piezoelectric properties) to 5.6%. Relative error was 1.1% for aluminium plates and 0.9% for polystyrene plates when considering error sign. Thus, it was proved that the proposed methodology for adjusting g31 variable as function of the vibration frequency and the stress level increased the accuracy of the analytical calculations. It is worthy to note that the implemented procedure was exactly the same than for the parameters used in the numerical simulations, so this procedure has been generalised for both calculation tools.
Finally, it was observed that the accuracy of the numerical simulations and the analytical calculations was in the same range. This fact proved that easier analytical calculation tools are powerful enough to deal with the modelling of cantilever piezoelectric devices.

Extrapolation for a theoretical case: analytical approach and literature comparison
Although the proposed procedure for calculating piezoelectric parameters from mechanical stress level and  Plot of r/r N vs. k 4 e to fit k 1 (slope) and k 2 (independent term) values of Eq. 1. Polystyrene specimens mechanical excitation frequency is limited to two substrate materials and low vibration frequencies in this research, its practical implications have to be highlighted to better understand its significance and to promote using it.
The proposed analytical methodology was applied to a theoretical case of energy harvester to evaluate the practical consequences of considering or not considering the stress and frequency influence on the values of the piezoelectric parameters for the calculation of electric energy production.
The considered theoretical case was a piezoelectric generator constituted by 10 bimorph cantilever devices cyclically activated by an imposed displacement whose amplitude was constant (1 mm) but excitation frequency decreased along time following an exponential curve. This generator was part of a complex energy harvesting device under development, whose conceptual design is shown in Fig. 6a. Every bimorph cantilever was constituted by a 93 mm 9 45 mm 9 0.5 mm aluminium (E = 68.9 GPa) plate. P-876.A12 piezoelectric patches were bonded on each side at 12 mm from the fixed plate ending. The fixed area covered a length of 10 mm. A 2 mm gap between the border of the fixed area and the piezoelectric patch were kept free for durability reasons. The imposed displacement was applied on a 3 mm length surface at the opposite edge of the plate. This configuration set a free cantilever length of 80 mm. These boundary conditions are represented in Fig. 6b.
Analytical and numerical approaches reached similar output accuracy in the previous analysis. Hence, only Experimental frequency and voltage reference, initial analytical output voltage, interpolated and adjusted piezoelectric constant and results of the adjusted analytical procedure and the corresponding relative error respect to experimental results Parameters used for the fitting of the piezoelectric properties in the analytical calculations analytical calculation was used for the theoretical evaluation of the energy harvester because of the lower complexity. Imposed displacement was characterised by the theoretical curve presented in Fig. 7, which started with a frequency of f 0 ¼ 8 Hz and finished when the frequency reached 1 Hz after 10 s. The frequency reduction followed the exponential law: f ¼ f 0 e À0:208t where t is time in seconds.

À3
Vm N ln Hz ð Þ ; k g 3 ¼ 3:1E À 3 Vm N ) was used to estimate and to automatically update the value of the piezoelectric voltage coefficient (g31) along the analytical calculation. An analytical calculation without considering the frequency and the stress influence and using piezoelectric voltage coefficient, g31 = 11.3E-3 Vm/N, provided by the manufacturer, was also carried out to compare with. Calculations considered bimorph configuration, which influenced beam stiffness, but the voltage of only one piezoelectric patch was calculated. The energy production of the full generator would be 20 times the one presented below.
Results of analytical calculations of both considered situations (manufacturer's g31 value and time-adaptative g a 31 value depending on frequency change and stress level) are presented in Fig. 8. It was observed that the voltage output amplitude was constant when a constant value of g31 coefficient was used (dashed line in Fig. 8a). In this case, the cumulative energy generated by one piezoelectric patch increased linearly with time (dashed line in Fig. 8b). In contrast, the results that considered the adaptative calculation of the g a 31 coefficient associated to stress and frequency dependency showed a decreasing voltage amplitude along time (solid line in Fig. 8a) resulting in a logarithmic like cumulative energy evolution (solid line in Fig. 8b).
For the considered reference time (10 s) of this theoretical case, cumulative energy production considering adaptative piezoelectric voltage coefficient reached 1/3 of the energy production estimated with the constant value of this variable provided by piezoelectric patch manufacturer.
To complete the analysis, it would be interesting to compare the proposed model with existing evidences. However, there are not publications, as per authors knowledge, that combine the analysis of different stress levels and different excitation frequencies below 10 Hz and provide the required data for applying the proposed

Conclusions
After performing experimental tests, analytical calculations and numerical simulations on cantilever bending piezoelectric devices with imposed cyclic displacement applied at frequencies below 10 Hz, the following conclusions were obtained: • Experimental output voltage amplitude increased when increasing the mechanical excitation frequency with fixed displacement amplitude and fixed deformation shape. Thus, properties that control mechanical-electrical relationship depended on vibration frequency. • Calculations considering constant piezoelectric parameters did not properly represent experimental evidences and brought relative errors up to 81% for numerical simulations and 392% for analytical calculations. • The newly proposed formulation to calculate piezoelectric parameters, which considered linear dependency on mechanical stress level and logarithmic dependency of imposed displacement frequency, accurately represented experimental results with a relative error around 2% for numerical simulations and 1% for the analytical approach. • This novel proposed formulation meets previously published evidences and it has been successfully generalised for two calculation tools: numerical simulations (adjusting parameters e31, ep33) and analytical calculations (adjusting g31). • In the case of low frequency vibration devices, considering the dependency of piezoelectric properties on the mechanical excitation frequency and the mechanical stress leads to significantly lower energy production estimation in comparison with the results obtained using the constant parameters suggested by piezoelectric manufacturers.
Finally, it is proposed to encourage researchers to perform tests at different stress levels and different low-range excitation frequencies so to collect evidences of the influence of these parameters on piezoelectric properties. the results of the experimental campaign. First author has written the article and proposed analytical methodology. Second author performed data post-processing, collaborated in discussion and revised the manuscript.
Funding Open Access funding provided thanks to the CRUE-CSIC agreement with Springer Nature. Partial financial support was received from COMSA company. Partial financial support was received from SORIGUÉ , S.A. (research contract PIEZOROAD).
Data availability Full data is available under request. Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http://creativecommons. org/licenses/by/4.0/.