A stiffness compensated piezoelectric energy harvester for low-frequency excitation

In this work, a stiffness compensated piezoelectric vibration energy harvester is modelled and tested for low-frequency excitations and large input amplitudes. Attracting magnets are used to introduce a negative stiffness that counteracts the stiffness of the piezoelectric beam. This results into a nearly statically balanced condition and makes the harvester a nonresonant device. A distributed parameter model based on modal analysis is used to model the output of the energy harvester. This model is extended by including the negative stiffness, endstop mechanics and force-displacement data to the model. The peak RMS power amounts 1.20 mW at 9 Hz and 3 g input acceleration. These are large inputs and serve to illustrate the case of having inputs larger than the device length. Furthermore, to benchmark the energy harvester in this work, the efficiency is evaluated in terms of generator figure of merit and is compared to prior art. This peak efficiency amounts to 0.567%, which is relatively large for its range of excitation. From the output that has been obtained with this design, it can be concluded that stiffness compensation can make a piezoelectric energy harvester competitive in terms of generator figure of merit at low-frequency excitation with input amplitudes exceeding the device length.


Introduction
Remote sensor networks are increasingly becoming important to monitor all sorts of processes [1]. Nowadays, batteries are commonly used to power them. However, they deplete over time and must be replaced, which can be challenging when the sensor is placed somewhere hazardous or hard to reach. As vibrations are omnipresent, transducing the vibrational * Author to whom any correspondence should be addressed.
Original Content from this work may be used under the terms of the Creative Commons Attribution 4.0 licence. Any further distribution of this work must maintain attribution to the author(s) and the title of the work, journal citation and DOI. power into useful electric power can provide a sustainable alternative [2].
In many cases a dominant frequency is present in the vibration signal, allowing a tuned resonating mechanism to be a good solution to extract power from the signal [2]. However, this only works well at higher frequencies. At low frequencies (below 10 Hz), the driving motion amplitude of a sinusoidal input rapidly increases for constant acceleration. As a resonator relies on amplifying the input amplitude, the length of the vibration energy harvester (VEH) then rapidly grows to dimensions unsuitable for its intended applications. So at lowfrequency excitation with large amplitudes, resonance quickly becomes impractical.
Several nonlinear techniques exist to possibly solve this. A tubular electromagnetic VEH was made with an increased electric damping to allow larger inputs from human motion [3]. A rolling mass in a circular cavity was used to circumvent the displacement limit of the proof mass [4]. Hard stops are also of use to limit the proof mass amplitude though it compromises the output power [5]. The principle of frequencyupconversion (FupC) utilizes a low-frequency oscillator that induces an impulse response in a high-frequency oscillator that generates the output [6,7]. In [8] it was shown that superharmonic resonances can be used for low-frequency harvesting, enabling downscaling of devices. Bistable VEHs have shown to provide decent output power figures in combination to wide bandwidths [9][10][11]. However, the nonuniqueness of the large-power orbits and strong dependency on initial conditions form a barrier to reliably generating a large output [12]. Stiffness compensation can also be used for low-frequency energy harvesting and allows for lowering the resonance frequency [13,14]. In the limit, a zero stiffness device can be created [15].
However, though these methods are promising and all have their (dis)advantages, apparently it seems hard to unify three certain properties into one energy harvester: a small dimension in the excitation direction, operation at low frequency and high efficiency. Often the largest dimension of the VEH is in the excitation direction [3,16,17]. This is beneficial to output power as it scales quadratically with the internal displacement limit and linearly with the other dimensions [18]. Furthermore it was shown that the output power has a cubic dependency on excitation frequency, so a small length VEH excited at low frequency tends to perform poor. For electromagnetic transduction, the detrimental effect of downscaling also influences this [19]. Furthermore, when a high efficiency, defined by the generator figure of merit (FoM g ) [20] is to be found at a low frequency it requires a greater length of the VEH as in [16,21,22]. Last, when a small length in the excitation direction is sought in conjunction with high efficiency, the device tends to operate at higher frequencies [23].
The previous discussion is characterized by three keywords: 'small in excitation direction', 'efficient' and 'low frequency'. Figure 1 shows that these do not coexist that easily. Two characteristics can be picked, and their intersection results in the opposite of the characteristic that is left. The intersection of all three, a VEH small in the excitation direction operating efficiently at a low frequency with large inputs seems infeasible up till now.
The research objective of this work is to end up in the intersection of all three characteristics by utilising stiffness compensation. A small piezoelectric cantilever beam is used which is generally too stiff to have resonance at low frequencies [24]. Therefore, the approach in this work is to compensate most of the stiffness of the piezo through addition of negative stiffness. This brings the piezo close to static balance and creates a nonresonant device. The negative stiffness can be formed by bistable flexure mechanisms [25,26]. Others have used magnets in attracting or repelling configuration [10,27]. Attracting magnets are used in this work. Due to the lowered stiffness, large deformations can be obtained more easily. This is beneficial for power generation and efficiency. A piezoelectricbeam claims little space in the excitation direction, allowing a more planar design. Combining this with the expected increase in efficiency at low frequency it can be a good candidate for ending up in the intersection of the three characteristics. The dynamics and output voltage of this prototype will be modelled and experimentally validated. Its performance will be compared to prior art in literature.
The next sections are organized as follows: in section 2, two important parameters are discussed that are used to benchmark the harvester performance with respect to prior art. The design of the harvester is discussed along with modelling approach. In section 3, the results from the simulation are compared to experimental results. These are discussed in section 4 along with the performance of the harvester with respect to literature. Finally, the main conclusions are summarized in section 5.

Motion ratio and generator figure of merit
Two parameters that are of major importance throughout this work are the motion ratio λ and the generator figure of merit FoM g introduced in [20]. These are defined as: Where L z is the length of the VEH along the excitation direction, Y 0 is the input motion amplitude, P RMS is the RMS output power, ρ pm is the density of the proof mass, V is the package volume of the VEH and ω is the radial frequency corresponding to the driving motion. The motion ratio is the ratio between the length of the VEH and the applied driving motion, see figure 2. In the ideal case where the frame claims no space, this length would be the displacement limit of the proof mass. It is a measure to determine to what extent the VEH can be used as a resonator. For instance, it can be used as a resonator when λ > 1 as in this case, resonant amplification of the input motion is possible. Hence, a resonator with a large Q-factor has a large motion ratio as well. When λ ⩽ 1, the input motion is larger than the length of the VEH and resonant amplification is not possible, resulting into a nonresonant device.
The generator figure of merit FoM g is an improved version of the volume figure of merit FoM v discussed [18,28]. It uses the density of the proof mass and instead of using V 4/3 , it uses VL z as the power output is more dependent on L z than the other directions [18]. This principle is shown in figure 2. The upper tubular configuration with larger L z like [3] is often seen in literature due to its larger potential for power generation. The lower configuration with smaller L z is more planar and is focussed on in this work for obtaining a better efficiency at a low motion ratio. The FoM g enables a fair performance comparison between both configurations by removing the excitation length bias.

Mechanical design
To compensate the stiffness of the piezoelectric beam, attracting magnets have been chosen in the configuration shown in figure 3. Details on the specifications can be found in table 1. By using attracting magnets in this setup, the introduced negative stiffness can be tuned through adjusting the distance between the magnets. As can be seen from figure 3, the frame acts as an endstop to limit the motion of the proof mass. A few reasons for this are: (1) the strain needs to be limited as to not damage the piezo, (2) it prevents the proof mass magnet from clinging to the fixed magnets and (3) it is a mechanism to transfer momentum into the proof mass motion.
The VEH is designed to function around excitations below 10 Hz with input amplitudes that are larger than the device length meaning a motion ratio lower than one. Normally, to ensure that maximum power is extracted from the energy harvester, the load resistance is assessed by using the impedance matching criterion as in [29]. As was seen in [24], the natural frequency of a piezoelectric beam increases with the load resistance connected to it, which implies that the load resistance affects the stiffness. Therefore, the load is kept constant at 1 MΩ.
In order to compensate the stiffness of the piezoelectric beam, its stiffness needs to be known. Therefore, a forcedisplacement measurement was taken with the setup discussed in section 2.5 to assess the stiffness at the uncompensated condition with a load of 1 MΩ connected. This measurement can be seen in figure 8(a). Then, a COMSOL Multiphysics electromagnetic simulation was carried out to find the distances d m1 and d m2 between the proof mass magnet and fixed magnets to provide an appropriate negative stiffness. The proof mass magnet was displaced inbetween the fixed magnets by a parametric sweep and the force was evaluated, resulting in a force-displacement relation. A parametric sweep was applied to d m1 , d m2 and the proof mass displacement to find the negative stiffness that fits best to the positive stiffness of the beam. The stiffness has been verified with the setup in figure 6. The results are shown in figure 8(a).

Fabrication
The piezoelectric beam used in the energy harvester is an in series connected Morgan Ceramics PZT-508 bimorph of 46 × 6 × 0.76 mm. An N35 Neodymium block magnet has been fixed to its tip with epoxy glue to form a proof mass. The frame of the VEH is 3D printed with 100% infill in PLA.
As the stiffness compensation becomes increasingly delicate as the total stiffness approaches zero and to account for printing tolerances, the frame was printed a few tens of millimeters larger in the L z direction. It was then manually sanded down to obtain the correct magnet distances.

Model
To model the dynamics and voltage output of the stiffness compensated piezo, a distributed parameter model is preferred over a lumped parameter model as an accurate description of the strain distribution is necessary. Furthermore, it was shown that a distributed parameter model could be used to provide corrections to a lumped model [24,30]. The main question is what this strain distribution looks like. One could reason that if no stiffness is reduced, modal analysis is suitable if resonance occurs. On the other hand, if the stiffness is fully removed, it could be reasoned that the deflection pattern can be described by the static deflection pattern of a cantilever beam with a force and moment at its tip multiplied by some temporal forcing function: Where δ is the tip deflection, P the tip load, L the beam length, x the beam length coordinate, YI the bending stiffness, M the tip moment and f the temporal forcing function. The eigenfunction from modal analysis and the static deflection pattern form the bounds where the real deflection pattern should be in.
To compare them, they are plotted in figure 4 along with their second derivatives which are directly related to the bending strain. All functions have been normalized. The figure shows that the deflections are very similar and that the strain diverges with the beam length. However, near the clamping the largest strain can be found and here the error is minimal. It can therefore be assumed that the deflection pattern is representative and modal analysis can be a representative method for modelling. Furthermore, modal analysis is a well-developed method and therefore convenient to use.
The modelling approach is as follows. First, modal analysis is used for the case of a cantilever beam. Next, the endstops are included in the equations and then the negative stiffness is added. Finally, the force-displacement data of the compensated beam is added to the equations. This last step makes sure that the remanent stiffness and hysteresis are included in the model. These steps will be discussed next.

Modal analysis.
In modal analysis, the relative displacement of the beam is formed by a summation of vibration mode shapes obtained at the resonance frequencies. In this case, only the first vibration mode is used as it closely describes the deflection pattern obtained in this harvester. The comprehensive distributed parameter model of Erturk and Inman is used for the model [24,31]. In their work, they had found closed-form solutions for the output voltage. In this work, nonlinearities are included through the addition of endstops and hysteresis. Therefore, nonlinear techniques based on energy methods can be used as an alternative [10,[32][33][34]. However, due to the numeric nature of the hysteresis measurement included in the equations of motion, the coupled equations in modal coordinates for an in series connected piezo are used in this work and solved numerically in MATLAB. Their derivation is concisely shown. The nomenclature is kept the same for easy reference. In figure 5, a schematic diagram of the piezo is shown. To transform the deflection w to modal coordinates the following infinite series is used: Where ϕ r (x) is the mass-normalized eigenfunction corresponding to eigenmode r and η r (t) is the temporal response corresponding to that same mode. The eigenfunction is defined as: where ς r = sin λ r − sinh λ r + λ r Mt mL (cos λ r − cosh λ r ) cos λ r + cosh λ r − λ r Mt mL (sin λ r − sinh λ r ) .
The modal eigenvalue λ r is found by solving the characteristic equation: 1 + cos λ r cosh λ r + λ r M t mL cos λ r sinh λ r − sin λ r cosh λ r − λ 3 r I t mL 3 cosh λ r sin λ r + sinh λ r cos λ r + λ 4 r M t I t m 2 L 4 1 − cos λ r cosh λ r = 0.
Where M t , m, L and I t are the tip mass, beam mass per length, clamped beam length and tip inertia, respectively. To find the amplitude C r , the eigenfunction ϕ r (x) is normalized by the following orthogonality condition: The bending stiffness is defined as: Where b, Y s , c E 11 , h p and h s are the beam width, Young's moduli and layer thicknesses of the substrate and piezo layer, respectively. Using the bending stiffness, the eigenfrequency can be found as: The coupled beam and electric circuit equations confined to the first mode in modal coordinates are as follows: where the forcing term F 1 (t) and electromechanical coupling θ are equal to: Where ζ 1 , V, C p , R l , e 31 , Y 0 and ω are the damping ratio, output voltage, single layer capacitance, load resistance, piezoelectric constant, driving motion amplitude and driving motion radial frequency, respectively. Most parameters can be found in table 1. The details on the derivation of these equations can be found in [24,31].

2.4.2.
Endstops. The endstops are defined as an additional spring and damper placed in parallel [5,35,36]. Figure 5 shows a schematic interpretation of the endstop. The force encountered by endstop contact can be described as: Here F s is the endstop force, K s the endstop stiffness, w s is the amplitude after which the stopper is hit and C s is the endstop damping. The endstop stiffness K s is found by Hertz contact mechanics. As the Hertz contact stiffness behaves nonlinearly with indentation, its stiffness value was raised until its corresponding indentation matched the indentation found in the final simulations. The endstop damping was assessed experimentally by a shaker test at 16 Hz and 2 g. By measuring the proof mass displacement and velocity with laser sensors, the coefficient of restitution was found, which can be translated to the endstop damping [36].

Negative stiffness.
To add negative stiffness, the negative stiffness from the finite element simulation or the measurement could be used. However, as indicated, the following step is to include the hysteretic force-displacement measurement at the compensated state. The introduced negative stiffness never exactly equals the stiffness of the piezo, leaving a remanent stiffness. This remanent stiffness is also present in this force-displacement measurement. As a result, this remanent stiffness would then be present twice in the model. Therefore, it is chosen to numerically add a negative stiffness that fully compensates the positive stiffness through subtraction of a tip deflection force: By fully removing the stiffness and using the remanent stiffness present in the force-displacement measurement of the compensated beam, an accurate description of the stiffness in the compensated state can be obtained.

Hysteresis.
As the stiffness has been fully erased in the model to provide a 'clean slate', the remanent bending stiffness and the hysteresis are to be found. These are found by quasistatically measuring the force-displacement of the compensated beam with the setup in figure 6. The measurement is taken for a connected load resistance of 1 MΩ. In figure 8(b) the measured hysteresis is shown in a force-displacement diagram and is named F H . The data has been filtered with a Savitzky-Golay filter and has been fitted with a cubic interpolant in order to implement it into the model. The implementation of the force-displacement data is as follows: at peak displacement, one of the endstops is hit. This hit is detected as an event in the MATLAB ODE-solver. After the hit, the hysteretic force switches from curve, e.g. the left endstop is approached from the red curve in figure 8(b) and the force switches to the blue curve afterwards and stays on it until the right endstop is hit. This is representative as long as peak-peak motion is present. A smaller loop is shown as well in figure 8(b) in black and green. This is for the a specific case when no peak-peak motion is obtained and the proof mass bounces from one endstop only. This case will also be validated.

Integration in the modal equations.
As the descriptions of the endstops, the negative stiffness and the forcedisplacement have been found, they can be integrated into the coupled modal beam equation (11). The resulting equation is found by equation (17). Note that the added terms are premultiplied by ϕ 1 (L) 2 . This is because the terms have to be converted to modal coordinates given by equation (4) and are premultiplied by the mode shape through mode shape orthogonality. The equations of motion are thus formed by equations (12), (13) and (17).
There is one thing that must be noted by taking this approach. As pointed out in Erturk and Inman [37], in a clamped-free beam, no tip force may be present at the free tip as this is one of the boundary conditions. In this model, the first vibration mode shape of the clamped-free condition is used and the added negative stiffness acting as a tip force theoretically violates the boundary condition. However, for now it is assumed that the mode shape will not be affected that much and that it may be negligible for the output, more about this in section 4.

Force-displacement setup
For the design and modelling of the VEH, force-displacement measurements are necessary. The setup to measure these is shown in figure 6. The setup consists of a Moons' 24Q-3AG stepper motor connected to an Almotion LT50-TR-G8-200 linear stage. A Futek LSB200 sensor is used to measure the force followed by a Scaime CPJ analog transmitter. The displacement is measured by a Micro-Epsilon optoNCDT1300 laser distance sensor and ELC DR07 resistance decade boxes are used to set the load resistance. The outputs are measured by an NI 9215 DAQ of which the data is analysed in MATLAB.

Harvester excitation setup
In order to validate the model, a custom air bearing stage with a linear motor is used to excite the energy harvester at low frequencies and large amplitudes. The setup that is used is shown in figure 7. An NI cDAQ-9174 with NI-9263 and NI-9215 modules is used. The NI-9263 generates the input signals that are sent to the stage which is feedback controlled. An ME-Meßsysteme AS28e accelerometer is used to check the applied acceleration and corrections to the input signal were made if the acceleration was not accurate enough due to positioning errors from the feedback loop. A resistive divider is used to lower the voltage from the piezo and is sent along with the output from the accelerometer to the NI-9215 to record the data.

Measurement procedure
In order to assess how accurately the simulation can predict the performance of the VEH, the RMS power frequency response is measured. In order to do so, a sinusoidal excitation is applied at an acceleration between 1.5 and 3 g and frequency between 2 and 10 Hz. Large accelerations are used to make sure the input is larger than the device length in most cases. The measurement is not taken as a sweep: at every single frequency, the measurement is started from zero initial conditions and the measurement time is adjusted per frequency as to reach steady-state conditions.

Benchmarking with respect to prior art
In order to benchmark the performance of the VEH in this work, it is compared to prior art found in literature. This is done by plotting the FoM g -λ space. This approach was also taken in [20], where a comprehensive comparison of VEH performance for different classes of VEHs was carried out. It must be noted that literature falls short on reporting all the variables needed to calculate the motion ratio λ and FoM g . In some cases however, the necessary parameters could be derived from other parameters or be estimated. For instance, in the case of a resonating cantilever, the size of the proof mass and the mechanical Q-factor are used to find L z .

Results
In figure 8 the force-displacement measurements are shown. Figure 8(a) shows the force-displacement of the uncompensated beam with a 1 MΩ load connected along with the simulated and validated negative stiffness. Figure 8(b) shows the force-displacement measurements for the compensated state of the VEH. A large loop is shown that indicates full range motion and a smaller loop as well for smaller motions. In figure 9 the simulated displacement, velocity and voltage along with the measured voltage is shown at the condition of 4 Hz and 1.5/3 g. Figure 10(a) shows the steady-state RMS power frequency response. This is done at a constant sinusoidal acceleration between 1.5 and 3 g and between 2 and 10 Hz. The simulation is compared to the measurements in terms of RMS output power at 2.5 and 3 g. Figure 10(b) shows the trend in output voltage and velocity. In figure 11, the performance of the energy harvester is evaluated in terms of motion ratio λ and FoM g and is compared to prior art from literature. Here, two classes can be seen, single degree-of-freedom (SDoF) and FupC.

Time domain
In figures 9(a) and (b), the blue line indicates the simulated displacement and velocity for the case where both endstops  are hit, which is the case for all measurements at 2.5 and 3 g. In this case, the measurements match the simulations best, as the used hysteresis loop is also based on full range motion. Figures 9(a) and (b) also show what happens if no full range motion is obtained, e.g. at 1.5 g. From the displacement in figure 9(a) it is seen that the proof mass detaches from the lower endstop and does not reach the upper endstop. Figure 9(c) shows that for the full range motion case, the voltage from the measurement matches well with the simulation. For the case of 1.5 g, shown in figure 9(d), extra measures must be taken. Here it can be seen that taking the full range motion hysteresis loop greatly overestimates the output voltage. To improve the simulation performance for this specific case, a hysteresis loop should be used that is representative for the lower displacement. This measured loop is shown in figure 8(b). Using this smaller loop as input for the simulation, it is seen that the voltage is resembled more closely. The consequences of this are however severe: for every separate non-full range motion, a separate hysteresis loop is necessary in order to predict the output within acceptable accuracy. This results into a large set of input data, but there is a method to reduce the dataset. In [38] it was demonstrated for piezoelectric actuators that inner hysteresis loops could be predicted by fitting curves to the outer loops and projecting them to the start and endpoints of the inner loop (black-green in figure 8(b)). So the small loop in figure 8(b) can be predicted by only knowing the large outer loop. This can be a good starting point to model non-full range motion, which will occur for lower input amplitudes or non-sinusoidal input motions.

Frequency response
The RMS power frequency response is shown in figure 10(a). The frequency response was measured at 1.5, 2, 2.5 and 3 g yet only for 2.5 and 3 g it is compared to the simulations as only in these cases full range motion is obtained. The peak performance is found at 9 Hz and 3 g where an RMS power of 1.20 mW is delivered. The effect of increased input acceleration is very clear. At 1.5 g only one endstop is hit, resulting in a low linear response. This makes sense as the output voltage depends on the strain rate integrated over the beam length, which is lower in this case [24]. At 2 g, the second endstop is hit slightly at some frequencies and at others it is not, resulting in a larger scatter in power. At 2.5 and 3 g both endstops are hit. These curves show an interesting behaviour. First, the response is linear until 6 Hz. This is also seen in other nonresonant devices such as in [39]. After 6 Hz, the response continues in a sublinear trend and a power saturation seems to appear.
The simulated output power response at 2.5 and 3 g is also shown in figure 10(a). It can be seen that at frequencies lower than 6 Hz, the simulated power closely predicts the measured. The error remains below 10%. However, from 7 Hz and onwards, the simulation rapidly diverges from the measurement and overestimates. The simulated voltage starts to overestimate, resulting in a rapid increase in RMS power.
The exact reason behind the rapid error growth is not fully understood yet. The results suggest that the exclusion of the negative stiffness in the mode shape function ϕ r (x) is admissible; the difference in power is small at frequencies lower than 6 Hz. If the mode shape function were not representative, a larger error should be seen for all frequencies, as in the simulation the deflection pattern remains constant for all frequencies.
Nevertheless, including the magnetic forces in the boundary condition could improve the performance.
The growth of the error in power can also be seen from the output voltage. Figure 10(b) shows the peak voltage in the measurement and simulation and the simulated velocity. It can be seen that the measured voltage stops to increase after 6 Hz, whereas the simulated voltage monotonically increases along with the simulated velocity. This is important, as the velocity directly relates to output voltage to the deformation rate as shown in equation (12). It is therefore most likely that at the higher frequencies energy is lost resulting in a lower velocity. This could be due to additional damping terms such as strain rate damping [40] or nonlinear losses in the endstop collisions. Furthermore, it is seen in piezoelectric actuators that hysteresis loops between input voltage and displacement show a rate dependence [41,42]. If this principle also translates between force and displacement, the used force-displacement loop measured at quasistatic condition could become less representative at the higher frequencies and additional terms may be necessary for compensation. Further development by including such terms is subject for future work.

Performance comparison to prior art
In figure 11, the measured performance of the energy harvester is compared to prior art in terms of FoM g and λ. From the plot it can be noted that prior art is inclined to larger motion ratios. The FupC harvesters tend to be less efficient than SDoF. The VEH of this work has a relatively good FoM g performance compared to other works within its range. The FoM g peaks at 0.567% at a motion ratio of 0.18. Although this is still low, larger values can be expected as no optimization has been applied yet. Nevertheless, the data show that having a design that uses a piezoelectric beam in combination with stiffness compensation results into a competitive efficiency at low frequencies.
Several things can be noted from the performance belonging to this work. First, the FoM g is relatively constant Figure 11. Generator figure of merit (FoMg) performance comparison of the energy harvester in this work to prior art as a function of motion ratio [3, 4, 6, 7, 10, 16, 22-24, 39, 40, 43-52]. SDoF stands for single degree-of-freedom, FupC for frequency-upconverter. throughout λ at constant acceleration. For instance, at 2.5 g a factor of 25 can be found in λ but only a factor of 1.38 in FoM g . To put this in perspective, resonating VEHs without endstops were analysed in [20] and showed a factor 16 in λ and 344 in FoM g . This clearly demonstrates the operational differences between resonant and nonresonant devices.
From figure 11, it is seen that the peak FoM g does not correspond to the peak power or peak acceleration. The peak FoM g is found at a motion ratio of 0.18 which corresponds to a frequency of 3 Hz and acceleration of 2 g. So although a higher frequency or input acceleration results into more output power, this clearly does not benefit the FoM g efficiency. This implies that after a certain value for the acceleration or input frequency, although the output power increases, the FoM g does not. Hence, there is a certain threshold after which the VEH starts making less practical use of its physical dimensions and properties.

Conclusion
In this work, the stiffness of a piezoelectric VEH is compensated through addition of attracting magnets. It was conceived in order to provide an energy harvester that has a small length in the excitation direction, that operates at a low frequency with large inputs and that is more efficient. A distributed parameter model from literature based on modal analysis was further extended by the addition of endstops, negative stiffness and hysteretic force-displacement data. For peakpeak displacements, the RMS power difference between simulation and measurement has been analysed. Between 2 and 6 Hz, the error between simulation and measurement remained below 10%, yet it rapidly grows when the frequency is further increased until 10 Hz. Therefore, a modal analysis based distributed parameter model can be found to be promising for modelling the dynamics and output of a cantilever-based stiffness compensated piezoelectric beam. However, work is necessary to improve the simulation performance over a wider range of frequencies.
The measured RMS peak power was obtained at 9 Hz and 3 g and was equal to 1.20 mW. When it comes to performance in terms of generator figure of merit (FoM g ), a peak value of 0.567% was obtained at 3 Hz and 2 g. The excitation levels are high but serve to illustrate the case of input motions larger than the device length. It can be concluded that stiffness compensation of piezoelectric beams can enable energy harvester designs to have a competitive efficiency at low-frequency excitation with large input amplitudes.

Data availability statement
All data that support the findings of this study are included within the article (and any supplementary files).