Optimization of piezoelectric cantilever energy harvesters including non-linear effects

This paper proposes a versatile non-linear model for predicting piezoelectric energy harvester performance. The presented model includes (i) material non-linearity, for both substrate and piezoelectric layers, and (ii) geometric non-linearity incorporated by assuming inextensibility and accurately representing beam curvature. The addition of a sub-model, which utilizes the transfer matrix method to predict eigenfrequencies and eigenvectors for segmented beams, allows for accurate optimization of piezoelectric layer coverage. A validation of the overall theoretical model is performed through experimental testing on both uniform and non-uniform samples manufactured in-house. For the harvester composition used in this work, the magnitude of material non-linearity exhibited by the piezoelectric layer is 35 times greater than that of the substrate layer. It is also observed that material non-linearity, responsible for reductions in resonant frequency with increases in base acceleration, is dominant over geometric non-linearity for standard piezoelectric harvesting devices. Finally, over the tested range, energy loss due to damping is found to increase in a quasi-linear fashion with base acceleration. During an optimization study on piezoelectric layer coverage, results from the developed model were compared with those from a linear model. Unbiased comparisons between harvesters were realized by using devices with identical natural frequencies—created by adjusting the device substrate thickness. Results from three studies, each with a different assumption on mechanical damping variations, are presented. Findings showed that, depending on damping variation, a non-linear model is essential for such optimization studies with each model predicting vastly differing optimum configurations.


Introduction
Over the last decade there has been a growing increase in research in the field of vibrational energy harvesting. The principle is to convert ambient vibrational energy into electrical energy with practical devices primarily targeted at replacing or supplementing the batteries used to power wireless sensors for condition monitoring, e.g. a tyre pressure sensor [1]. The use of piezoelectric materials is seen to be the most effective transduction method [2], reflected by the fact of piezoelectric devices receiving the most interest. For general information on vibrational energy harvesting the reader is referred to a compendium in the form of a recent book [3]. However, more specifically for piezoelectric energy harvesters, the reader is directed to Piezoelectric Energy Harvesting [4].
There has been a great deal of work on the linear modelling of cantilever piezoelectric harvesters. Sodano et al [5] used the Rayleigh-Ritz procedure to estimate the power output from a cantilever mounted piezoelectric generator. Ertuk and Inman [6] later developed a distributed-parameter electromechanical model for energy harvesters. Validation was also completed by the authors, with theoretical results in good agreement with the experimental data [7]. Patel et al expanded on this by accurately incorporating the effects of non-uniform beams created by altering piezoelectric layer coverage. Utilisation of the model showed that vast improvements in performance are achievable, with the experimental data indicating the modelʼs validity [8]. However, despite the level of interest in the area, reliable nonlinear modelling of piezoelectric harvesters validated through detailed experimental work is generally lacking.
Models exist to simulate energy harvesters operating in conditions or arrangements where non-linearity is induced intentionally. This approach is utilised to assist in overcoming a major limitation to vibrational energy harvester usability-a limited operational bandwidth [9]. In early research, Ramlan et al [10] theoretically investigated the effects of introducing non-linearity via spring hardening with system behaviour represented in the form of the Duffing equation. System bandwidth was found to increase, with the peak magnitude of the generated voltage remaining on a similar level to that predicted by a linear model. Similar theoretical and experimental research, whereby non-linearity is induced through the use of magnets, can be found in [11,12], to name a few examples. More recent advances on the use of magnets are presented by Tang et al in [13,14]. The authors propose the use of a magnetic oscillator in place of a fixed magnet with experimental results indicating a 100% improvement in bandwidth and a 42% improvement in peak power at an acceleration level of 2 m s −2 over both fixed magnet and linear harvester designs [13]. The piezoelectric energy harvester-magnetic oscillator design is realised through the inclusion of an additional cantilever with a magnetic proof mass, above the conventional energy harvester cantilever. Readers are directed to [14] for detailed information, through experimental work, on the use of magnets in improving the functionality of energy harvesters through monostable/bistable device designs as well as designs targeted at frequency upconversion for use in low-frequency applications. An alternative method for intentionally inducing non-linearity is through the use of a static axial pre-load. Masana and Daqaq [15] accurately modelled and experimentally validated such a scenario for piezoelectric energy harvesters. Findings indicated that axial pre-load can be beneficial in the improvement of both device bandwidth, by introducing non-linearity, and peak power, by increasing the electrical damping.
Non-linearity naturally exists in vibrating systems undergoing large deflections (geometric non-linearity), and when certain materials are involved in harvester composition, particularly associated with the piezoelectric layer. As shown by Joshi in [16], non-linear constitutive equations can be used to define material non-linear behaviour, with Crespo da Silva and Glynn [17] providing fundamental work on geometric non-linearity in beams. Few works can be found which attempt to incorporate these effects into the modelling of energy harvesters, e.g. [18][19][20], with experimental validation providing confidence in the developed models. Stanton et al provide the ground work for the modelling of inherent piezoelectric material non-linearities in energy harvesters [19]. Model validation is also provided through experimental testing on an off-the-shelf bimorph device. The model is restricted to the analysis of uniform samples, with symmetry through the thickness, where harvesters are comprised of layers identical in length. The effects of substrate material non-linearity and geometric non-linearity are not considered in [19]. The representation of such factors is seen to be important in the modelling of certain types of energy harvester, namely thin film flexible energy harvesters. Mutsuda et al [21] have proposed an ocean power generator which comprises underwater hanging highly flexible harvesters attached to an elastic floating unit. The floating unit can generate electric power from wave oscillations and wave breaking whilst hanging units utilise vibrations created by underwater currents, vortices and oscillation. Naturally, as a result of low oscillatory frequencies occurring in an ocean environment (≈0.8-1.2 Hz), the harvester must be manufactured from flexible materials. Silicone rubber for the substrate layer(s) and polyvinylidene flouride (PVDF) for the piezoelectric layer(s) are the materials of choice. The usage of such flexible devices is not limited to the ocean environment and can be utilised in, for example, wind flutter applications, extracting energy on/near bridges. The materials used in the manufacture of such devices, and the large deflections which they undertake, indicate the requirement for an improved analytical non-linear model. Our study is heavily motivated by these devices. It is important to be able to predict how nonlinear behaviour will affect the frequency response of the harvester (the extent of resonance shift and peak magnitude reduction) in order to ensure the harvester is designed to operate most efficiently at the dominant excitation frequency.
In this paper a versatile non-linear model will be developed which has the ability to accurately predict piezoelectric energy harvester performance. The model will include a transfer matrix sub-model to determine the dynamic behaviour of segmented beams which are generated when altering the coverage of the piezoelectric layer. Coupling with energy expressions for segmented harvester structures will allow for the optimization of devices in terms of their piezoelectric layer coverage. This was not considered before, e.g. in [18] or [19], but is seen to be essential when designing performance/ cost efficient devices. Substrate material non-linearity will also be included in the development of the model. This factor is highly important as the model is expected to be used on thin film flexible devices in the near future. Additionally, the curvature of the beam required to incorporate geometric nonlinearity will be obtained through differential geometry. The paper is laid out in the following manner. Section 2 will disclose information on the modelling approach and the embodiment of material and geometric non-linearities. Hamiltonʼs extended principle along with the calculus of variations will be used to obtain the final two equations of the motion (transverse displacement and voltage). This will be followed by a substantial section on model validation through experimental work, section 3. Results obtained from four different samples manufactured 'in-house', with various lengths of piezoelectric material will be presented here. Section 4 illustrates the importance of using the non-linear model during device optimization over a linear model in certain operating conditions. Finally, the paper will end with closing remarks and suggestions for future work.

Modelling a piezoelectric harvester system including various non-linearities
In this section the details behind the theoretical modelling of energy harvesters are presented, including both geometric and material non-linearities. To increase model versatility by allowing for alternations to piezoelectric layer coverage, a component of the model will utilise the transfer matrix method [22]. The accurate prediction of natural frequencies and mode shapes will result from such an approach. Electrical aspects of the system will consist solely of a load resistor, with earlier work on more realistic circuitry, for example, storing generated energy in a capacitor, provided previously in [8].

The energy harvester
The general methodology behind the model development is similar to that undertaken by several other researchers, i.e. [1,23]. Material non-linearity is introduced in the form of additional terms found in constitutive equations, and geometric non-linearity will result from an inextensible beam assumption. The extended Hamilton principle, along with the calculus of variations, is used to obtain equations of motion in the time domain.
A schematic of the structure is provided in figure 1. The harvester considered consists of a composite two-layer Euler-Bernoulli beam, with piezoelectric material perfectly bonded to a substrate layer. x 1 represents the distance of the piezoelectric layer from the clamped end, and x 2 is the length of the piezoelectric layer. The current Euler-Bernoulli beam assumption, in which shear deformation and rotary effects are neglected, is reasonable for the vast majority of expected harvester geometries, and justified by Dietl et al in [24]. Other notation in figure 1 is taken in conjunction with [23], with the Newtonian inertial co-ordinate system represented by x y z ( , , ) and the local co-ordinate system represented by ξ θ ζ ( ) , , .
2.1.1. Constitutive equations. The constitutive equation of a material is used to relate axial stress, σ 11 , to axial strain, ε 11 , and in addition, for the piezoelectric material, the electric displacement, D 3 , to axial strain. Herein, index 3 refers to the y-direction, i.e. through the thickness of the material, and index 1 refers to the x-direction, i.e. parallel to the beam length. For the piezoelectric material [23] one has  where superscript p refers to the piezoelectric layer, E p is the piezoelectric material Youngʼs modulus, and E field is the electric field strength. d 31 is a piezoelectric material constant, and ε 33 is the material permittivity. The difference between the above equations (equations (1) and (2)), and the linear relationships used in [7] for example, is the inclusion of higher order terms with constants μ 1 and μ 2 . These terms represent non-linearity, and are both specific and unique to each 'batch' of piezoelectric material.
For the substrate material one has where superscript s refers to the substrate layer and E s is the substrate material Youngʼs modulus. Note how a coefficient of material non-linearity, μ s1 , has also been included for the substrate layer. As a result of experimental setup uncertainties, the importance of this term is realised in section 3.1. However, as alluded to in the introduction, its inclusion is paramount for the modelling and performance predictions of certain harvester compositions, work which is not presented here.
Equations governing system response will be determined through application of the Hamilton extended principle. This requires both the Lagrangian of the system and work done on the system, where the Lagrangian refers to the difference between kinetic energy, T, and potential energy, U. An expression for each of these terms is obtained in the following subsections.   where L is the length of the substrate layer, or + + x x x 1 2 3 , s d is the small element length, and A s and A p are cross-sectional areas for the substrate and piezoelectric material, respectively. Terms v s t ( , ) and u s t ( , ) refer to the transverse and longitudinal deflection respectively, see figure 1, while ′ ( ) donates the derivative with respect to the arc length, s. Note how, herein, for the deflection terms, the independent variables s and t are excluded for ease of reading. The expression for EA(s) must take into account non-uniform material distribution, realised through the use of Heaviside functions: s s s p p p 1 1 2 Strain in the beam can be expressed in terms of the distance from the neutral axis, y, and beam curvature, ρ: 11 where ρ is obtained using differential geometry, see for example [25], and given by: with higher powers approximated as: 3 3 Utilising equations (1)-(3) and equations (6)- (8), in addition to the assumption that the electric field is uniform throughout the piezoelectric thickness, , transforms equation (4) to: where the terms K 1 through K 5 are given by:  where y is the location of the neutral axis from the bottom of the substrate layer. I s , I s1 , I p1 and I p2 are provided by: , where t s is the substrate thickness, and b s and b p refer to the substrate and piezoelectric layer widths, respectively. The inextensibility condition is used to eliminate independent longitudinal vibrations, u, in equation (10), by relating them to transverse vibrations, v. For the inextensibility condition to be satisfied, the strain along the neutral axis [23], ϵ 0 , must equate to zero, where ϵ 0 is given by: Expansion, rearrangement and utilisation of Taylorʼs expansion leads to the following relationship: 2 Using the above relationship in equation (10) results in a two-degree-of-freedom system, one dependent co-ordinate for displacement and another for voltage. The next step is to obtain the system Lagrangian which requires an expression for the kinetic energy.
2.1.3. Kinetic energy and the system Lagrangian. The kinetic energy of such a system can be represented as: Note that the transverse-longitudinal relationship has been used to eliminateu. m(s) is the mass per unit length and can be defined as: where in turn ρ s and ρ p are the substrate and piezoelectric material density, respectively. The Lagrangian of the system can be defined as usual: which upon substitution of equations (10) and (15) leads to: In the following subsection, the extended Hamilton principle is utilised in order to obtain the two equations of motion.

Governing equations of motion.
Before the calculus of variations can be used, an expression for external work done on the system, , requires formulation, and is a combination of the base excitation,ẅ t ( ) b , and the electric potential energy: is the electric charge generated by the energy harvester. Applying the extended Hamilton principle, i.e. and utilising the calculus of variations, leads to the following equations of motion:  with the following associated boundary conditions: A classical modal analysis technique, namely the Bubnov-Galerkin method, is used to obtain simplified ordinary differential equations from the existing partial differential equations previously shown. Using this approach, the beam deflection, v s t ( , ), is expressed as an infinite sum of the products of normalised eigenvectors, W s ( ) r (see section 2.2), and time-dependent generalised co-ordinates, η t ( ) where 'r' refers to the mode number. Substituting this relationship into equations (21) and (22), in addition to using orthonormality conditions, yields the following two  and Note how proportional damping, γ r , has been introduced in equation (25), to accommodate for energy dissipation from the system. ω r refers to the natural frequency of the r th mode and is obtained using the transfer matrix method (see section 2.2). In obtaining equation equation (26), the rate of change in charge,q t ( ), i.e. current, has been expressed as −V t R ( ) load . This formulation is acceptable as these works will simulate and present results from energy harvesters connected directly to a load resistance. C n r 1 to C n r 14 are the resulting constants independent of time, provided in the appendix.
Equations (25) and (26) represent the behaviour of unimorph vibrational piezoelectric energy harvesters under non-linearity inducing conditions. They can be solved simultaneously to determine transverse vibrations along the structure and voltage generated by the energy harvester. Numerical solving of the equations is achieved using ODE solvers in MATLAB ® , through the SIMULINK ® interface [26], with data recorded once a steady-state response has been achieved.

Transfer matrix model for segmented structures
The model presented in section 2.1 requires knowledge of the natural frequencies and mode shapes of the harvester. In this section the transfer matrix method [22] is used to obtain this information for a segmented cantilever beam, by taking into account the length and position of the piezoelectric layer. As shown in figure 2 the beam is split into three sections. The substrate material alone makes up sections 2 and 3, while section 2 comprises both piezoelectric and substrate materials.
The notation used to define the mechanical forces and deformations at the modes of each element are shown in figures 2 and 3.
The exact beam function for the transverse motion of the i th section of a segmented beam is given by: where a, b, c and d are constants which are ordinarily determined using boundary conditions, x is the distance from the left side of the beam segment and l i is the beam segment length. β i can be defined by: where E i is the element Youngʼs modulus, I i is the element area moment of inertia, ρ i is the element mass density and A i is the element cross-sectional area. The middle section is a composite piezoelectric/substrate beam and the properties of this element can be calculated using the equivalent flexural rigidity: s ys p yp composite where the parallel axis theorem can be used to obtain I ys and I yp , such that:   The mechanical forces and deformations on the righthand side of any beam element, z Ri , can be related to the lefthand side, z Li , through the use of a transfer matrix, ( ) l U i i [22], where,  3 0 Constants C 0 to C 3 are defined by: Full details on how the transfer matrix and related constants are derived is provided in [22]. For the general case, where piezoelectric material (of a length shorter than the substrate layer) is centrally located on the beam, the overall transfer matrix of the system, U overall , is obtained from: can be used to obtain a matrix whose determinant yields the natural frequencies of the system. The 2 × 2 matrix of interest is extracted from the bottom right corner of the U overall matrix. Values of ω which produce determinant values of zero provide the natural frequencies of the system. Since mode shapes can be arbitrarily scaled, once natural frequencies are known the corresponding mode shapes are readily obtained by assuming one of the variables, i.e. clamped end shear force F, to be unity. Following this the mode shapes are scaled as required, i.e. to the mass of the structure, allowing for the validity of orthonormality conditions used in section 2.1.4.

Validation through experimental testing
In this section validation of the theoretical model presented in section 2 is provided through the experimental testing of uniform and non-uniform samples. The energy harvesters, manufactured 'in-house' and of unimorph type, are comprised of an aluminum substrate layer and a lead zirconate titanate (PZT) piezoelectric layer. Adhesion is realised through a combination of DP460 epoxy and a small amount of silver conductive epoxy, with complete details behind the consistent manufacturing procedure found in [27]. In terms of the testing procedure, samples are mounted in a clamp attached to a Data Physics GW-V4 electromagnetic shaker providing base excitation. A Stanford Research Systems SR785 dynamic signal analyser is used to output a harmonic signal to the shaker via a standard amplifier. In addition to this, the analyser has the capability to record two unique signal inputs. One channel was always used to monitor base acceleration, in order to keep it fixed during frequency sweeping, while the second channel provided measurement for either voltage across a resistor or tip velocity. The base acceleration was measured using a PCB Piezotronic accelerometer, model number 352C23, and the tip velocity was measured with a PolyTec OFV-055 laser vibrometer, with velocity readings readily converted to displacement.
The coefficients of material non-linearity, i.e. μ 1 and μ 2 , are unique to each 'batch' of piezoelectric material, and are not provided on the data sheets of material manufacturers. Such coefficients, along with μ s1 , are estimated using curve fitting techniques, as is standard practice in the field. Note that since all testing is performed close to fundamental frequencies during theoretical simulation, the inclusion of solely the first vibrational mode, in equations (25) and (26), is necessary, i.e. = r 1. The fixed dimensions, i.e. layer widths and thickness, along with relevant material properties are provided in table 1. Each sample in this section will be subjected to four increasing levels of base excitation-0.5 ms −2 , 2.5 ms −2 , − 5 ms 2 and 7.5 m s −2 .The validation of the model is provided by estimating the coefficients from one sample alone and using the obtained magnitudes in the theoretical model to predict and compare the results for other samples. Sample 1 (and a sample of only Al) was used to obtain μ 1 , μ 2 and μ s1 . The obtained magnitudes were then used to predict the behaviour of non-uniform samples 1 and 2. The predicated behaviour will be shown to be in good agreement with experimental trends, indicating that the model is valid for a wide range of piezoelectric material coverages.

A sole substrate layer
Firstly, the testing of a sole substrate layer, i.e. a device comprising of no piezoelectric material, was undertaken in the experimental setup. The results from this test are provided in figure 4, superimposed with data from theoretical simulations. The sample used had an extended length of ± 43 1 mm, where 'extended length' refers to the overhang length from the clamp face.
The experimental data in figure 4 indicates a shift in the resonance from 295.9 Hz to 294.1 Hz when subjecting the sample to a base acceleration of 0.5 ms −2 and 7.5 ms −2 , respectively. The magnitude of this shift can be accommodated by setting μ s1 to − × 2 10 13 Pa; found through curve fitting in MATLAB ® [26]. The importance of including an additional form of non-linearity in conjunction with geometric and piezoelectric material non-linearity is evident from the example presented here. The softening phenomena would not have been reproducible had substrate material non-linearity been excluded from theoretical model formulation. The observed non-linearity could have resulted from other sources, such as boundary condition non-linearity or inertial nonlinearity (see [28] for a detailed review on non-linearities in vibrating systems). However, the inclusion of substrate material non-linearity, albeit in a crude manner, will be shown in section 3.2 to provide excellent experimental-theoretical agreement across a range of harvester samples. As has been observed by previous researchers [29], the magnitude of mechanical damping experienced by the structure is dependent on the magnitude of acceleration, and found to be approximately linear across the tested range. Note that this is not in agreement with the non-linear damping assumption made by Stanton et al in [19]. Damping magnitudes corresponding to theoretical plots in figure 4 were extracted from experimental data using the half-power points method [31] and are shown in figure 5.
In depth model development incorporating the variation in mechanical damping is beyond the scope of this work.
However, for the experimental setup utilised, one can briefly include air flow damping, stick-slip at the clamped end and material damping, as the most likely causes of the observed variations. Further details on variation in mechanical damping can be found in [30].

Uniform and non-uniform energy harvesters in closed circuit conditions
Following on, in this section the outcome of testing three energy harvester samples is presented, each comprising differing piezoelectric layer lengths. The substrate and piezoelectric layer lengths for each sample are provided in table 2. Note that for the non-uniform samples, the piezoelectric material length is reduced from the free end, i.e. clamping occurs on both piezoelectric and substrate layers, or = x 0 1 in the theoretical model.
The results in this section are from three different samples tested in closed circuit electrical conditions, i.e. = R 0.
load This eliminates one of the piezoelectric coefficients of non-linearity, namely μ 2 , facilitating the curve fitting process. Note that the magnitude of the substrate material coefficient of non-linearity, μ s1 previously found in section 3.1, is used herein when generating theoretical results. Experimental data for the uniform conventional sample, superimposed with theoretical results, is provided in figure 6. Figure 6 shows that substantial non-linearity does exist in real situations and a softening phenomenon is witnessed once again. This observation will be a continual theme indicating that material nonlinearity dominates over geometric non-linearity for this type of cantilever piezoelectric energy harvester. Numerically, the maximum deflection equates to 0.0017% of the overhang length, which reinforces the reasoning behind why geometric non-linearity is negligible for this particular device. The resonant frequency of the structure was found to be 360.4 Hz at is 0.48%. Note how, although peak magnitude is well predicted by the model, the off-resonant behaviour is not. This is possibly due to imperfections in bonding during the sample manufacture, with results from samples with shorter piezoelectric layers (figures 8 and 9) showing better off-resonance matching. The variation in mechanical damping is again seen to resemble a quasi-linear increase with base acceleration, with data provided in figure 7 for completeness.
Results from samples comprising shorter piezoelectric layers will now be shown in order to demonstrate the versatility of the model. Figures 8 and 9 show results from the experimental testing and theoretical simulation of non-uniform samples 1 and 2, respectively. It is important to note that no more curve fitting for non-linearity terms is undertaken here. The theoretical data in each case is obtained using the , respectively. Figures 8 and 9 show that good experimental-theoretical agreement is also obtainable for non-uniform samples using previously determined magnitudes for material non-linearity coefficients. The shift in resonant frequency between the two excitation extremes in each case is found to be 4 Hz and 1.6 Hz, increasing with piezoelectric material length. Prediction of the resonant frequency by the theoretical model, at an acceleration of 7.5 m s −2 , results in percentage errors of 0.06% and 0.29%. This highlights the accuracy and versatility of the model used in section 2 in predicting the behaviour of piezoelectric harvesters subjected to excitation conditions that induce non-linearity. The mechanical damping was again found to increase in a quasi-linear manner with results presented in figure 7. The occurrence of differing rates of increase in damping ratio with base acceleration (26%, 70% and 94% for the uniform sample and non-uniform samples 1 and 2, respectively) emphasises the difficulties in predicting damping magnitudes. In the following section, the energy   harvester will be connected to a load resistor in order to determine the magnitude of the remaining non-linear coefficient, μ 2 .

Energy harvester connected to an electrical load
The magnitudes of μ 1 and μ s1 found from previous closed circuit testing are still valid as the sample tested here was manufactured from the same 'batch' of sheet piezoelectric material. The two variables which require estimation are μ 2 , found in equation 2, and γ r = r ( 1), found in equation 25. Dimensions for the sample under consideration are provided in table 1 with a ± 45.02 mm 1 mm extended length. The magnitude of mechanical damping (0.0062, 0.0085, 0.0108 and 0.0125), corresponding to the results in figure 10 at the four excitation levels, was obtained through closed circuit testing and curve fitting as previously demonstrated. Figure 10 provides the results obtained when a 150 kΩ resistor is introduced to the system. The first observation to note is that good agreement between experimental and theoretical voltage responses is achievable without need for the non-linear coefficient μ 2 . Following this observation, the effects of μ 2 on the theoretical results are investigated, with the realisation that μ 2 has little influence on theoretical frequency responses. This is another difference when comparing the proposed theoretical model with that developed by Stanton et al in [19], where the electro-elastic non-linear constant is utilised. It is believed that this coefficient can be assumed to be zero in the majority of energy harvesting scenarios due to the inherently low voltage levels. Applications utilising piezoelectric material for actuation are subjected to higher voltage levels and in these situations the non-linear coefficient,   μ 2 , would have a more significant impact on the theoretical results. As a final note, in terms of energy harvesting, the vision for the devices being proposed by Mutsuda et al [21] is to generate power on the kilowatt scale and so analysis on these devices would eventually require the predictions and utilisation of μ 2 . From comparing the frequency shift for cases with and without a resistor, as one would expect, the level of experienced non-linearity reduces when energy is extracted from the system by the addition of an electrical load. The frequency shift (from an excitation of 0.5 m s −2 to 7.5 m s −2 ) reduces to 3.5 Hz from a 6.5 Hz shift obtained when the load resistor is excluded. Around 390 Hz, there is noise in the experimental data in figure 10. This is due to a resonant frequency from the clamp and shaker arrangement interfering with vibrations of the energy harvesting device.
The non-linear constants obtained in this paper are only valid for this particular 'batch' of substrate and piezoelectric material. For each new 'batch' of material, experimental testing must be undertaken to determine the magnitudes of the three material non-linearity coefficients. It is advised that the following be undertaken to achieve this.  • It is advisable that testing of only substrate cantilevers be performed initially. This will eliminate all but one material non-linear coefficient, μ s1 , facilitating the fitting process. Tip displacement FRFs should be used here. • Following this, it is suggested that a complete energy harvester sample (either uniform or non-uniform) be tested in closed circuit conditions. The fitting process should now be performed to obtain μ 1 , using the already determined μ s1 magnitude. Again tip displacement FRFs should be used here. The fitting process can also be used here to obtain variations in mechanical damping with base acceleration. • Finally, the testing of a complete energy harvester connected to the load resistor can be undertaken. Voltage FRFs can be used in this case, with curve fitting used to provide the user with the only remaining non-linear coefficient, μ 2 .

Linear and non-linear model comparisons during device optimization
In this section, device optimization is performed in relation to coverage of the piezoelectric layer. Theoretical results from both a linear and non-linear model will be provided and comparisons made. The width of both layers is taken to be 5 mm and the thicknesses of the substrate and piezoelectric layer are taken to be 0.67 mm and 0.5 mm, respectively. The material properties assigned to the device are given in table 1.
The length of the substrate layer will remain fixed at 50 mm with the piezoelectric layer length varying from 50 mm to 1 mm. Length reductions are made from the free end, i.e. x 1 remains 0 in figure 1. Note how the Euler-Bernoulli assumptions made during the modelling procedure hold true as the overall beam length remains at 50 mm and only the length of the piezoelectric material is being reduced.

Creating devices with identical fundamental frequency
Before performing the study on piezoelectric coverage, a discussion on why a constant fundamental frequency approach is used is essential. Altering one geometric parameter alone is not advisable as this creates devices with differing fundamental frequencies. Suppose that the conventional device (both layers being the same length) is designed for an application to operate most effectively at F Hz, the dominant excitation frequency. Any changes to the piezoelectric layer length, whilst keeping all other geometric parameters constant, will alter the fundamental frequency creating an ineffective device. To avoid this mismatch between fundamental and dominant excitation frequencies across all designs, in this work the thickness of the substrate layer is used as a control parameter. Through this simple procedure, unbiased design comparisons can be made. Figure 11 shows the substrate thickness required for each piezoelectric layer length to create configurations with identical fundamental frequencies.
The data plotted in figure 11 is obtained by sweeping through a range of substrate thickness for each design, plotting fundamental frequency vs thickness, and interpolating to find an accurate value. The trend is as one would expect. As material is reduced from the free end, the mass of the device reduces at a greater rate than its stiffness. This causes an increase in fundamental frequency which must be countered by reducing the substrate layer thickness. After a critical point, which depends on the materials/initial size/inclusion of a tip mass, etc., the stiffness begins to reduce at a higher rate than the device mass and so increases in the substrate thickness are required. For the non-linear model, non-linear constants, μ s1 , μ 1 and μ 2 obtained in section 3.2 are used during the study.
The effect of variations in mechanical damping during optimization studies is highly important and cannot be ignored. Detailed work on changes in mechanical damping, and predictions during piezoelectric coverage optimization using a linear theoretical model, can be found in [27]. For the purpose of this work three scenarios will be considered.
• The damping magnitude will remain constant across all configurations. • Damping data from both the uniform and the second nonuniform sample (PZT coverage of 100% and 21% respectively) is used to obtain a power relationship between the device volume and damping ratio. Therefore in this scenario it is assumed that the magnitude of mechanical damping depends on the device volume where reducing the volume by approximately 25% reduces the damping by 50%; see the dashed line in figure 12 for numerical values. • Experiments in [27] showed that predicting trends between mechanical damping and the PZT length is extremely difficult. The author found that, overlying the general trend, large variations in damping occurred between samples, variations which were predominantly due to mounting and inconsistent clamping force. One observed trend between the PZT length and mechanical  Figure 11. Plot indicating the substrate thickness required for designs with various piezoelectric layer lengths. All configurations have identical natural frequencies. Figure 12. Two possible trends between the mechanical damping ratio and piezoelectric material coverage (note that the substrate thickness is also changing in accordance with figure 11). The first trend, represented by the dashed line, assumes a power relationship between damping and volume. The second, represented by the solid line, is obtained by scaling the findings found in [27].   damping showed a linear increase in damping with length to approximately 66% coverage, thereafter plateauing off for longer PZT lengths. By assuming the major contributor to the damping magnitude to be the amount of adhesive, a similar trend can be used in this work by scaling according to data from the uniform sample in figure 7. The final trend between the damping ratio and the length of PZT which will be used in the optimization study is shown in figure 12. Figure 13 clearly shows differences in the performance trend obtained whilst using either a linear or non-linear model when a constant damping assumption is made. The linear model suggests that maximum voltage is generated when the piezoelectric material covers the entire beam, whereas the non-linear model suggests maximum voltage is generated for a device with a ≈5 mm long piezoelectric layer. For cases where the piezoelectric layer is short the difference between the model outputs is minimal and a linear model would suffice. Due to the volume of piezoelectric material the effects from material non-linearity, which is the dominant effect, are reduced. However, as the material length increases, so too does the extent of non-linearity. This is the cause of the observed discrepancies for devices with longer piezoelectric coverage. Note how for devices with the thinnest substrate layers, i.e. piezoelectric lengths of approximately 30 mm from figure 11, large differences between the model predictions exist. Naturally, the deflections experienced by these devices will be greatly increased which induces non-linear behaviour, thereby resulting in the relatively low voltage levels seen in figure 13.
In figure 14 results from both linear and non-linear models, while assuming that damping varies through a power relationship with volume, are presented. The general trend is clearly very different from that obtained in figure 13 where damping was assumed constant. Peak performance is seen to occur for devices with PZT coverage of approximately 68% when utilising the linear model, and 100% coverage whilst utilising the non-linear model. For devices which exhibit low mechanical damping there is large divergence between the two models as a result of increased non-linearity in devices. In this case, devices with reduced PZT coverage are showing poor levels of performance due to the presence of large mechanical damping. Information on damping ratio magnitudes utilised here can be found in figure 12.
When the damping magnitude is assumed to partly vary with adhesive length, in the manner shown in figure 12, both models predict the same optimum configuration; see figure 15. In this scenario, the samples with increased PZT coverage experience higher mechanical damping suppressing motion, thereby reducing the effects of non-linearity on peak performance. For this reason the differences between results from a linear model and a non-linear model are much smaller than those obtained when assuming a low constant damping ratio, figure 13. The low mechanical damping for samples with reduced lengths of PZT is responsible for the observed spike in peak voltage. Similar trends, i.e. high performance from devices with short PZT layers, have previously been reported through experimental work in [27].
This brief study highlights the importance of using a nonlinear model in combination with detailed knowledge on mechanical damping variations during device optimization, if warranted by the excitation conditions. Although for the last tested case both models predicted the same optimum PZT coverage, the non-linear model is essential in predicting the resonant frequency at which peak performance occurs; see figure 16. The device performance would have suffered by 5% if this shift was unknown to the energy harvester designer. In this example the shift, and detriment to performance, is small because the optimum configuration only has a 10% piezoelectric coverage, however, this will not always be the case.

Conclusion
A robust non-linear model to predict the dynamic response of piezoelectric cantilever energy harvesters has been developed and presented. Since piezoelectric material is known to behave in a non-linear fashion, even at moderate excitation levels, the inclusion of piezoelectric material non-linearity is seen as essential. Material non-linearity was also included for the substrate layer in addition to geometric non-linearity for the beam. The inclusion of substrate material non-linearity was important in terms of extending the applicability of the model and accounting for effects observed during experimental testing. Geometric non-linearity is realised by applying an inextensibility condition, and material non-linearity is incorporated through the addition of higher order terms in constitutive equations. In order to ensure piezoelectric layer coverage can still be optimized, a transfer matrix model was developed allowing for accurate predictions of the eigenvalues and eigenvectors of segmented structures. Detailed model validation is provided through the use of harvesters manufactured 'in-house'. This is then followed by a comparison of results from the non-linear model with those of a linear model while predicting the performance of energy harvesters.
Closed circuit testing was initially undertaken, allowing for simplification of the fitting process when obtaining the magnitude of the non-linear material constants. For this 'batch' of piezoelectric material a μ 1 magnitude of −7 × 10 14 Pa was found to provide adequate matching between the experimental and theoretical data. This was seen to be the case across a range of samples, each comprising a different piezoelectric layer length. The devices were then attached to a load resistor in order to determine the magnitude of the second non-linear constant, μ 2 . It was concluded that the inclusion of this non-linear constant is not required when modelling harvester devices, possibly due to the low voltage levels associated with energy harvesting in comparison to using piezoelectric material for actuation.
To determine the magnitude of the material non-linearity coefficient for the substrate layer, tests were performed in the absence of piezoelectric material. The magnitude of μ s1 was found to be −2 × 10 13 Pa, which, as one would expect, is lower than that of the piezoelectric material. Another important finding was that mechanical damping increases in a quasi-linear manner with the base acceleration. This means that simple scaling, i.e. double the excitation to double the performance, does not strictly apply. Harvesters must be simulated in realistic operating conditions in order to accurately predict how they will perform.
A case study on device optimization has also been presented here. Three different assumptions were made in regards to the mechanical damping variations with PZT length: (i) a constant damping ratio for all devices, (ii) a power relationship between damping and device volume, and (iii) a relationship obtained from a similar experiment undertaken in [27], extrapolated for this scenario. The results indicate that trends between piezoelectric layer coverage and performance are highly dependent on damping variations. The use of a non-linear model for estimating the optimum configuration in terms of peak power is debatable due to uncertainty in the mechanical damping prediction. However, the model is essential for acquiring knowledge on the extent of changes in the resonant frequency for the optimum configuration.
In terms of future work, experimental testing on harvesters with other material compositions is recommended. This will assist in gauging the applicability and limitations of the developed non-linear model. Testing on highly flexible devices such as those designed to be used in ocean applications, proposed by Mutsuda et al [21], is currently in progress. Silicone rubber and PVDF materials make these devices more sensitive to geometric and substrate material non-linearities, which will reinforce the usefulness of the developed model.