Analytical and Experimental Investigation of a Curved Piezoelectric Energy Harvester

Piezoelectric energy harvesters have traditionally taken the form of base excited cantilevers. However, there is a growing body of research into the use of curved piezoelectric transducers for energy harvesting. The novel contribution of this paper is an analytical model of a piezoelectric energy harvesting curved beam based on the dynamic stiffness method (DSM) and its application to predict the measured output of a novel design of energy harvester that uses commercial curved transducers (THUNDER TH-7R). The DSM predictions are also verified against results from commercial finite element (FE) software. The validated results illustrate the resonance shift and shunt damping arising from the electrical effect. The magnitude, phase, Nyquist plots, and resonance frequency shift estimates from DSM and FE are all in satisfactory agreement. However, DSM has the advantage of having significantly fewer elements and is sufficiently accurate for commercial curved transducers used in applications where beam-like vibration is the predominant mode of vibration.


Introduction
Vibration energy harvesting involves the scavenging of ambient kinetic energy by transforming it to electrical energy using electromagnetic, electrostatic or piezoelectric devices [1][2][3][4][5]. Due to their effectiveness, light weight and small size, piezoelectric devices are extensively used as an alternative to batteries in low power devices [6]. Piezo materials can be integrated into lightweight structures, such as beams and plates, to convert electrical power into mechanical energy (actuator mode) or mechanical power into electrical energy (sensor mode or vibration energy harvesting mode). Considerable research into the modelling and analysis of piezoelectric devices has been performed with regard to actuators and sensors for smart vibration control, e.g., [7,8], vibration energy harvesting, e.g., [9][10][11][12][13][14][15][16], or even simultaneous vibration control and vibration energy harvesting, e.g., [17,18]. The focus of the present paper is vibration energy harvesting.
A review of vibration energy harvesting literature [19] shows a number of publications researching configurations designed to access the maximum voltage via various mechanisms, e.g., induced by moving loads on civil structure applications [12], vehicles travelling on a bridge [20], and aeroelastic galloping [21]. Sodano et al. [13] conducted an experimental comparison of the power generation capabilities of a variety of piezoelectric materials. With regard to mathematical modelling techniques, although numerical methods were used as early as 2002 [14], significant progress on modelling the electromechanical coupling was made from 2008 onwards following the analytical approach of Erturk and Inman [9,11,22,23].
Erturk and Inman [9] utilized the Euler-Bernoulli beam theory along with the constitutive equations of piezoelectric material to set up a distributed parameter model of a base-excited energy harvesting cantilever made up of a unimorph (i.e., single piezoelectric layer bonded to a metallic substrate). This model was then transformed using a modal transformation wherein the clamped-free modes with no electrical effect (short circuit conditions) were used to represent the flexural vibration of the beam relative to the vibrating base with a resistive shunt circuit. The resulting equation of motion (considering only one mode) was solved to obtain the displacement, voltage, and power response for harmonic excitation over a range of frequencies. The same authors later performed a similar analysis of base-excited cantilevered bimorph with tip mass, which was additionally validated by experimental investigation [9]. Both [9,22] illustrate the vibration damping effect and shift in response frequency induced by the electrical coupling.
Rafique and Bonello [24] also used the above-described analytical modal analysis method (AMAM) of distributed parameter energy harvested cantilevers to model a bimorph, this time without tip mass, and obtained satisfactory correlation with experimental results for voltage, displacement, and power frequency response functions (FRFs). For a more thorough validation, Nyquist plots of the FRFs were considered in addition to the more usual magnitude FRF plots. In [25], Dalzell and Bonello extended the AMAM approach to an experimentally validated analysis of an energy harvesting cantilevered bimorph shunted by an energy storage circuit comprising a diode in series with a capacitor for half-wave AC/DC rectification.
In another study, Bonello and Rafique [26] introduced the dynamic stiffness method (DSM) for the analysis of base-excited cantilevered unimorph and bimorph devices with or without tip mass. The model was based on the Euler-Bernoulli beam model with piezoelectric coupling and its results were verified against those obtained with the AMAM. The DSM can be extended beyond simple uniform-section cantilever beams more readily than AMAM. It was shown in [26] that the DSM could be applied to more complex systems composed of an assembly of straight piezoelectric beam segments with various boundary conditions using matrix assembly techniques similar to those used in the finite element (FE) method. However, unlike FE or AMAM, the DSM uses the exact function to describe the segment's vibrating shape under harmonic excitation. Hence, for the purpose of frequency domain analysis, the DSM offers an efficient computational approach in which the accuracy is independent of element size, unlike FE [27,28].
In most of the aforementioned research, the energy harvesting devices have taken the typical form of base-excited cantilevered straight beams with/without a tip mass. However, there is a growing body of research into the use of curved piezoelectric transducers for energy harvesting purposes. Such curved transducers were originally designed as actuators, and compared to their traditional counterparts, these curved devices not only offer more flexibility and reliability, but are also capable of being applied to harvest energy from lowfrequency vibrating structures [10]. Examples of curved piezoelectric transducers that have been researched are RAINBOW [29], LIPCA [30], and THUNDER [10,31], with the latter being considered the most promising due to its unique performance features compared to traditional unimorph, bimorph, and straight extensional actuators [32]. THUNDER (thin layer unimorph ferroelectric driver) consists of stainless-steel substrate, the piezo-electric wafer, and an aluminum top layer and adhesive bonding layers. Its curved shape is the result of prestress induced during the manufacturing process of the device. After processing in an autoclave, a curved prestressed device with high out-of-plane displacement is formed, in part due to a mismatch between the coefficients of thermal expansion corresponding to the PZT and the metal layer [32,33]. Mossi et al. [34] showed that THUNDER devices' performances can be significantly affected by design parameters, such as the geometry, number of layers, and the thickness of the layers. The same conclusion for other types of THUNDER were also observed in [35].
In 2016, Wang et al. [31] investigated how to maximize energy harvesting efficiency by changing the radius of curvature of a THUNDER as a tuning parameter. The first part of their study involved an experimental and analytical investigation for the output voltage and power across a load considering both standard AC and full wave AC/DC rectification. The thin shell device was modelled as a curved laminated beam and it was further assumed that the curvature was sufficiently low so that the curvilinear coordinate along the neutral surface could be used interchangeably with the rectangular coordinate. Moreover, the analysis represented the vibrating shape using a modal transformation approach based on only one mode (i.e., a Rayleigh-Ritz procedure). The effect of changing the radius of curvature of the THUNDER was investigated in the second part of [31] but no experimental validation was provided.
In 2017, Bharat Kathpalia et al. [10] studied, analytically and experimentally, the implementation of a commercial arched piezo energy harvester (precisely a THUNDER configuration) for smart paver tiles. Closed-form expressions for the voltage output to force input frequency response function (FRF) and power output to force input FRFs, for a simple resistive load, were derived and experimentally validated for range low-frequency (<10 Hz). It was concluded that tens of microwatts of power could be produced by a single arch at low force levels, and milliwatt levels of power can be achieved by using several arches.
In 2018, Hasan et al. [36] established an FE method to address the analysis of piezoelectric coupled fields. The main idea behind this research was to study the effect of changing the radius of curvature to see how it changes the output voltage and power of the THUN-DER device. They reported a notable agreement in experimental data and modal analysis carried out by ANSYS. A more recent study, carried out by Thonapalin et al. [37] in 2021, discussed the influence of temperature on different types of THUNDER samples, such as THUNDER 6R, 7R, and 8R. The experiments revealed, among other things, that high temperature badly affects the energy harvesting capability of the devices, particularly in the low frequency regime (e.g., reduction of 40% at 80°C) [37].
Motivated by the above-described advantages of the DSM, the novel contribution of the present paper is an analytical model of a piezoelectric energy harvesting curved beam based on the DSM and its application to predict the measured output of an energy harvester that uses commercial curved transducers (THUNDER TH-7R). The DSM predictions are also verified against results from commercial FE software. In the DSM approach the transducer is modelled as a curved laminated beam but, unlike [31] (where the Rayleigh-Ritz approach was used), no limitations on the degree of initial curvature are assumed. The analysis will include all layers of the transducer, including the adhesive layer, which was omitted in energy harvesting work [27,32] despite earlier evidence of its influence on actuation [29]. The energy harvesting device is designed such that the input excitation is applied to one end of the parallel transducers in the longitudinal direction, vibrating a tip mass attached at the other end. This design is novel for energy harvesting purposes, but a similar design has been used as an adaptive tuned vibration absorber by Bonello et al. [38] wherein the THUNDER devices were used in actuator mode to control the curvature and hence the tuned frequency.

Modelling by Dynamic Stiffness Method
The linear constitutive law is incorporated into the equations of motion of a laminated curved beam with constant radius of curvature. Assuming harmonic vibrations, the dynamic stiffness method (DSM) is then applied to obtain the dynamic stiffness matrix of a curved piezoelectric beam. The DSM method is applied to the energy harvesting device in Figure 1, which consists of two THUNDER devices pinned together at both ends with a tip mass, whose electrodes are connected in parallel. Utilizing a matrix assembly process enables the inclusion of attachments to the curved beams, and boundary conditions.
The material properties and the dimensions of the THUNDER device are given in Tables 1 and 2, respectively. The tip mass in Figure 1a is 410 g, including the mass of radial bearings at the free end. The material properties and the dimensions of the THUNDER device are given in Table 1 and Table 2, respectively. The tip mass in Figure 1a is 410 g, including the mass of radial bearings at the free end.

Equations of a Curved Piezoelectric Beam Segment
Although THUNDER is actually a curved thin plate rather than a curved thin beam, the way the transducers are set up in the energy harvester design (Figure 1), allows each to be modelled as a curved thin beam with a good approximation.
The strain for a curved thin beam is as follows.
where and are respectively the tangential and radial deformations at mid-surface, is the radius of curvature at mid-surface, is the curvilinear coordinate along the midsurface, and z is the distance between the evaluation point and the mid-surface in the thickness direction (see local coordinate axes in Figure 2). Stress and charge density can be expressed as follows [39].

Equations of a Curved Piezoelectric Beam Segment
Although THUNDER is actually a curved thin plate rather than a curved thin beam, the way the transducers are set up in the energy harvester design (Figure 1), allows each to be modelled as a curved thin beam with a good approximation.
The strain for a curved thin beam is as follows.
where u and w are respectively the tangential and radial deformations at mid-surface, r is the radius of curvature at mid-surface, s is the curvilinear coordinate along the mid-surface, and z is the distance between the evaluation point and the mid-surface in the thickness direction (see local coordinate axes in Figure 2). Stress and charge density can be expressed as follows [39].
where subscript k refers to layer no. k (see Table 1), k = p refers to the piezo layer (as per piezoelectric stress constant, and S 33 is permittivity at constant strain. To include material damping (c), the stress relations are modified as follows [26]: where subscript refers to layer no. (see Table 1), = refers to the piezo layer (as per Table 1 is the elastic modulus, ℎ is the thickness, is the voltage generated, is piezoelectric stress constant, and is permittivity at constant strain. To include material damping (c), the stress relations are modified as follows [26]: The equations of motions for a curved laminated beam in terms of displacements are given as: Note that represents the viscous damping coefficient and is mass per unit length. The procedure to derive these equations and the expressions for constants ( , , ) and operators ( , , and ) are given in Appendix A.
For Rayleigh-type damping the following relations are used in advance (assuming the material damping coefficient distribution through thickness follows the same profile as modulus of elasticity).

= ,
(6a) where and are defined as mass proportional and stiffness proportional damping coefficients.
Following the same procedure applied in [40], and can be related to modal damping ratios and and frequencies and of the first two modes of vibration, resulting in the relations. The equations of motions for a curved laminated beam in terms of displacements are given as: Note that c a represents the viscous damping coefficient and m is mass per unit length. The procedure to derive these equations and the expressions for constants (A, B, D) and operators (K, K, and χ) are given in Appendix A.
For Rayleigh-type damping the following relations are used in advance (assuming the material damping coefficient distribution through thickness follows the same profile as modulus of elasticity).
where α and β are defined as mass proportional and stiffness proportional damping coefficients. Following the same procedure applied in [40], α and β can be related to modal damping ratios ζ 1 and ζ 2 and frequencies ω 1 and ω 2 of the first two modes of vibration, resulting in the relations.
Since the experiments conducted (Section 4.2) cover only one mode (i.e., only ζ 1 and ω 1 are available) it is not possible to determine both coefficients α or β. In this case, one of the coefficients is arbitrarily set to zero and the other is fitted according to Equation (7a). In this paper α is set to zero (i.e., the stiffness proportional damping has been selected) and β is given by: Integrating the time rate of Equation (3) at midsection of piezoelectric layer (z = z pc ), and assuming that the electrode covers the extent of the segment (s = 0, s = S i ) leads to electrical current equation: where T 33 is permittivity at constant stress. Substituting Equation (1) into Equation (9a) yields the subsequent formula for the current: where Notice that Equation (11b) neglects the last term of Equation (1) since z/r for a thin beam is negligible.
Assuming harmonic variation in time of the generated voltage with frequency ω rad/s (i.e., AC voltage) and an external load of generic impedance Z(ω), Equation (10) can be transformed into the frequency domain using complex notation as follows where v is the complex amplitude of the harmonic voltage v which is then given by v = Re ve jωt where Re{} denotes the real part of {} and j = √ −1 denotes the imaginary unit.

Dynamic Stiffness Matrix of Curved Beam
Since the voltage is assumed to be harmonic, the vibration generating it is also harmonic, and the complex notation for the variables w and u in Equations (5) will be as follows u(s, t) = Re u(s)e jωt (13b) Substituting Equations (13) into Equations (5) and taking advantage of the separation of spatial and time variables in Equation (13), the homogenous part of the equations of motion will become as follows.
Now, assuming w = ae λs , u = γae λs and substituting in Equation (14), the eigenvalue problem is given by: The conditions for the existence of the nontrivial solution of Equation (15) requires that the determinant of the coefficient matrix of Equation (15) must be zero which yields a characteristic equation for wavenumber λ.
Solving Equation (16) for λ, the general solution of Equation (14) is shown as follows.
where (·) represents element-wise product operator. For a segment no. i, the solution must satisfy displacement and force/moment boundary conditions at both endpoints (s = 0 and s = S i , Figure 2). Definingd(s) = w u ψ T one can get, from Equation (17) and considering that ψ = ∂ w/∂s − u/r, the following expression.
Definingf(s) = V N − M T where N and M are the complex amplitudes of N and M (given by Equations (A1)) and V = ∂ M/∂s is the complex amplitude of the shear force V, one can substitute (17a) into (A1) to get In Equation (19): λ n = λ n 1 λ n 2 λ n 3 λ n 4 λ n 5 λ n 6 . v (i) is the complex amplitude of the voltage generated by the segment (no. i), derived by placing Equation (13) into Equation (12) and re-arranging: Substituting Equation (17) into Equation (20) v (i) = be 31 It should be noted that in the specific cases considered in this paper, the electrical load is a pure resistor R L , i.e., Z(ω) = R L . However, the use of the generic impedance term Z allows application to more complicated external circuits comprising resistors, capacitors, and inductances. Such complex circuits are useful in energy harvesting to maximize the power generated through electrical impedance matching [41,42], or even to achieve simultaneous energy harvesting and ambient vibration attenuation through a tuned shunted circuit [18].
With reference to Figure 2, the nodal degree-of-freedom vector d (i) and corresponding vector of forces and moments f (i) for segment no. i are defined as Finding a from Equation (22) and substituting it into Equation (23) where D (i) is the dynamic stiffness matrix of segment no. i. If lumped inertias of masses M i and M i+1 and moments of inertia I i and I i+1 are attached to end nodes nos. i and i + 1 of segment no. i, the following matrix is added to its dynamic stiffness matrix D (i) [26]: To avoid duplication of inertia, a given attached lumped inertia should be either added only to one segment, or shared (in any arbitrary proportion) between contiguous segments. Now considering T (i) to be transformation matrix between global and local coordinates ( Figure 2), where (θ i ) is the angle between tangent at node no. i and global x axis. The global dynamic stiffness matrix of the segment is therefore

Matrix Assembly and Frequency Response Functions
Due to symmetry of the energy harvester configuration in Figure 1a, the modelling need only consider one of the transducers, each carrying half the tip mass and shunted across twice the external circuit resistance R ext (i.e., in previous equations, Z(ω) = R L = 2R ext since transducers are connected in parallel across R ext ). With reference to Figure 1b, a total of seven segments were used to model the device. Five segments-2, 3, 4, 5, and 6-were used to model the region between the base joint and the tip joint (both joints coinciding with the centers of the radial bearings which act as pivots). Segments 1 and 2, and 6 and 7, correspond to the clamped regions and were therefore considered rigid and their density was chosen such that their moment of inertia (MoI) about their respective joint axis matches with that of the corresponding clamp, alternatively this clamp inertia could be applied using Equation (25) but for the sake of using a unified model for both DSM and ANSYS it is avoided; any slight excess clamp mass arising from choosing the density on the basis of its MoI is corrected by correspondingly adjusting the value of the tip mass ( Figure 1) which is included using Equation (25) with I i set to zero. Segments 3 and 5 correspond to regions of steel substrate extent that are not covered by the clamp-this gap region is necessary to avoid risk of cracking the piezoelectric material and is considered with edge length of 3 mm. Notice that the mid-surfaces of these short steel segments are not continuous with the mid-surface of the contiguous composite segment No. 4. This introduces a slight error in the model since the matrix assembly process described below assumes continuity of the mid-surface. However, the effect of this error is not considered significant (as indeed evidenced by verification against ANSYS in later sections) since the beam is thin. Using a matrix assembly process based on continuity of displacement at adjacent ends of contiguous segments, the dynamic stiffness matrix D of the whole system with 7 segments (Figure 1b) is obtained: (1) . . . (1) . . .
In Equation (28), d and f contain the complex amplitudes of the global nodal displacements at node nos. i (x i , y i , ψ i ) and the corresponding excitations The tip mass matrix (Equation (25)) is added to the rows and columns corresponding to node no. 7 (Figure 1b).
The boundary conditions are as follows: The receptance matrix α then relates d to f as follows where α is determined in the following way: • the rows relating to y 2 and y 7 , and columns relating to F y 2 and F y 7 are all padded with zeros due to the boundary conditions; • the remainder of α is α red = D −1 red where D red is obtained from D by eliminating the columns relating to y 2 and y 7 , and the rows relating to F y 2 and F y 7 . Now, with reference to Figure 1, the only effective non-zero excitation in f is F x 2 , since F y 2 are F y 7 are non-effective due to the boundary conditions (notice also that F x 7 = 0 since the inertia and damping effect from the tip mass at node No. 7 are considered in Equation (25) and the equivalent damping fit of Equation (8)).
Hence, the transmissibility of the device, relating the harmonic displacements/velocities/accelerations at tip (node No. 7) to the base (node No. 2) is given by: where α x 7 , F x 2 and α x 2 , F x 2 are the receptance matrix terms relating the displacements at nodes 7 and 2, respectively, with the transmitted force at the base. The voltage is produced by segment no. 4 (Figure 1b). Finding a from Equation (22) and substituting it in Equation (21) gives From Equation (20): where α d (4) , F x 2 is that sub-column of α relating the complex amplitudes of nodal degrees of freedom of segment no. 4, contained in d (4) , with the complex amplitude of the transmitted force at the base joint F x 2 F x 2 . Substituting Equation (33) into Equation (32), the frequency response function (FRF) relating v (4) to F x 2 is written as The voltage FRF relating v (4) to the base acceleration (normalised by g = 9.81m/s 2 ) is then written as The power FRF is calculated from voltage FRF using relations

Finite Element Modelling Using Ansys
Finite element (FE) modelling was performed using ANSYS (ANSYS Academic Research Mechanical, Release 19.1). The ANSYS elements that support the piezoelectric effect are PLANE223 (Figure 3a) and SOLID226 (Figure 3b). These were used in the piezoelectric layer. For the non-piezoelectric layers, they were replaced by their non-piezoelectric counterparts (PLANE183 and SOLID186, respectively). Both piezoelectric element types are compatible with CIRCU94 element, which is used to model external circuit (resistors in this case). In order to model electrodes of PZT patch, defined nodes on the upper surface and bottom surface of the piezoelectric layer are coupled into two separate sets. The PLANE223/183 elements lie in the x-y plane ( Figure 3a) and have unit depth in the z-direction. Therefore, the value of resistance should be multiplied by b (the width of the THUNDER), and the tip mass should be divided by b, while modelling with this (plane-type) element. Furthermore, for the model with PLANE223/183 elements, the Poisson ratio effect applies only in the x-y plane (the plane of the elements), i.e., does not apply along the z-direction. Hence, the model with PLANE223/183 elements shown in Figure 3a will be equivalent to a model based on beam elements in terms of deformation behavior (i.e., curvature in x-y plane only). type) element. Furthermore, for the model with PLANE223/183 elements, the Poisson ratio effect applies only in the x-y plane (the plane of the elements), i.e., does not apply along the z-direction. Hence, the model with PLANE223/183 elements shown in Figure 3a will be equivalent to a model based on beam elements in terms of deformation behavior (i.e., curvature in x-y plane only).
(a) (b) For the purpose of investigating mesh convergence, a purely mechanical setup (i.e., with the external resistance short-circuited) was considered, and a modal analysis was conducted to obtain and compare undamped natural frequencies of the device in Figure  1 using different element types. In addition to aforementioned PLANE223/183 and SOLID226/186, the element type SHELL281 (which does not support the piezoelectric effect) was also considered for the purely mechanical exercise (there is no piezoelectric shelltype element in ANSYS). All of these elements use quadratic shape functions to interpolate field variables. MASS21 was used to model the tip mass. For the 3D models (solid or shell elements), a symmetry condition was applied in the z-plane. For the free vibration analysis, the boundary conditions were applied as shown in Figure 3, i.e., zero displacement in z-direction at all nodes, zero displacement in y direction at tip and base joints, zero displacement in x-direction at base joint. The values for the first four natural frequencies of the device using the different types of ANSYS elements are presented in Table 3 (columns 2-6). The corresponding modal displacement shapes of the transducer section of the device modelled using solid elements are depicted in Figure 4. Convergence was achieved for an element edge length of 2 mm in circumferential direction, and for plane and solid models one element was enough through thickness of each layer. For shell and solid models, 3 mm of element edge length was enough in z direction. The number of elements used in each model for plane, shell, and solid elements were 350, 1211, and 4565, respectively. Figure 4 shows that the transducer section of the device vibrates approximately like a beam in the first four modes, with curvature being mainly in the x-y plane (particularly the first three modes). With reference to Table 3, the plane element results (which effectively describe beam-bending behavior of the THUNDER in the x-y plane, as explained above) are consistent with both the shell and plane results since the ~3 to ~5% shortfall in the natural frequencies can be attributed to the effective bending stiffness of each layer of the shell and solid element models being 1 (1 − ) ⁄ times that of a beambending model, where is the Poisson ratio of the layer [43]. In fact, neglecting the adhesive layer, assuming an average Poisson ratio of 0.3, and considering that the natural frequencies are approximately proportional to the square root of the bending stiffness, the percentage difference in natural frequencies between the plane element model and the  For the purpose of investigating mesh convergence, a purely mechanical setup (i.e., with the external resistance short-circuited) was considered, and a modal analysis was conducted to obtain and compare undamped natural frequencies of the device in Figure 1 using different element types. In addition to aforementioned PLANE223/183 and SOLID226/186, the element type SHELL281 (which does not support the piezoelectric effect) was also considered for the purely mechanical exercise (there is no piezoelectric shell-type element in ANSYS). All of these elements use quadratic shape functions to interpolate field variables. MASS21 was used to model the tip mass. For the 3D models (solid or shell elements), a symmetry condition was applied in the z-plane. For the free vibration analysis, the boundary conditions were applied as shown in Figure 3, i.e., zero displacement in z-direction at all nodes, zero displacement in y direction at tip and base joints, zero displacement in x-direction at base joint. The values for the first four natural frequencies of the device using the different types of ANSYS elements are presented in Table 3 (columns 2-6). The corresponding modal displacement shapes of the transducer section of the device modelled using solid elements are depicted in Figure 4. Convergence was achieved for an element edge length of 2 mm in circumferential direction, and for plane and solid models one element was enough through thickness of each layer. For shell and solid models, 3 mm of element edge length was enough in z direction. The number of elements used in each model for plane, shell, and solid elements were 350, 1211, and 4565, respectively. Figure 4 shows that the transducer section of the device vibrates approximately like a beam in the first four modes, with curvature being mainly in the x-y plane (particularly the first three modes). With reference to Table 3, the plane element results (which effectively describe beam-bending behavior of the THUNDER in the x-y plane, as explained above) are consistent with both the shell and plane results since the~3 to~5% shortfall in the natural frequencies can be attributed to the effective bending stiffness of each layer of the shell and solid element models being 1/ 1 − ν 2 times that of a beam-bending model, where ν is the Poisson ratio of the layer [43]. In fact, neglecting the adhesive layer, assuming an average Poisson ratio of 0.3, and considering that the natural frequencies are approximately proportional to the square root of the bending stiffness, the percentage difference in natural frequencies between the plane element model and the shell or solid element models is expected to be 100 6, which agrees with the percentage differences for plane vs. shell and plane vs. solid in Table 3.  The last two columns of Table 3 respectively show the DSM results for the undamped natural frequencies and the percentage differences relative to the ANSYS plane element results. It is clear that the ANSYS plane element results are much closer to the DSM results than to the ANSYS shell or solid element results. This is expected since the DSM model is directly based on beam-bending, and the ANSYS plane element model describes a similar beam-bending behavior, albeit modelled indirectly through plane elements (resulting in the small differences in the last column of Table 3). It should also be noted that, as seen from Section 2, DSM is a frequency domain method used to derive FRFs. Hence, unlike FE (ANSYS), the undamped natural frequencies in DSM are not found by solving an eigenvalue problem since the final relation is a system of equations (Equation (30)) relating the nodal displacements and forces in the frequency domain (rather than an equation of motion in the time domain involving discrete mass and stiffness matrices). However, the undamped natural frequencies can be determined by DSM from the frequency locations of the peaks of the transmissibility FRF (Equation (31)) under zero damping conditions. For zero electrical load and no material damping, these resonance peaks would be infinite, and their frequencies are more accurately located by noting the zero crossing points of the denominator of Equation (31) (i.e., the FRF , ). It should also be noted that, although the DSM undamped natural frequencies are determined from the FRF, the DSM is still regarded as an analytical method since the FRF terms are based on the analytical solution (Equation (17)) of the dynamic equations of the curved beam segment.

Experimental Testing
In this section, the experimental setup is presented (4.1) followed by a description of the estimation of the damping in the energy harvesting device (4.2). The last two columns of Table 3 respectively show the DSM results for the undamped natural frequencies and the percentage differences relative to the ANSYS plane element results. It is clear that the ANSYS plane element results are much closer to the DSM results than to the ANSYS shell or solid element results. This is expected since the DSM model is directly based on beam-bending, and the ANSYS plane element model describes a similar beam-bending behavior, albeit modelled indirectly through plane elements (resulting in the small differences in the last column of Table 3). It should also be noted that, as seen from Section 2, DSM is a frequency domain method used to derive FRFs. Hence, unlike FE (ANSYS), the undamped natural frequencies in DSM are not found by solving an eigenvalue problem since the final relation is a system of equations (Equation (30)) relating the nodal displacements and forces in the frequency domain (rather than an equation of motion in the time domain involving discrete mass and stiffness matrices). However, the undamped natural frequencies can be determined by DSM from the frequency locations of the peaks of the transmissibility FRF (Equation (31)) under zero damping conditions. For zero electrical load and no material damping, these resonance peaks would be infinite, and their frequencies are more accurately located by noting the zero crossing points of (ω)). It should also be noted that, although the DSM undamped natural frequencies are determined from the FRF, the DSM is still regarded as an analytical method since the FRF terms are based on the analytical solution (Equation (17)) of the dynamic equations of the curved beam segment.

Experimental Testing
In this section, the experimental setup is presented (Section 4.1) followed by a description of the estimation of the damping in the energy harvesting device (Section 4.2). Figure 5 shows an annotated photograph of the energy harvesting device and Figure 6. shows the experimental setup. The experimental setup was comprised of the following equipment.

•
A PC-controlled data acquisition (DAQ) system (not shown in Figure 6) consisting of hardware (LMS Scadas 5) and spectral analysis software (LMS Test.Lab). • Two accelerometers (PCB 352C22), one attached to the base of the energy harvesting device and the other attached to its tip mass. • A signal conditioner/amplifier, to remove the noise from signals received from accelerometers and amplify them before channeling to the DAQ system. • An electromagnetic (EM) shaker on which the energy harvesting device was mounted. • A power amplifier to control the gain of the random excitation signal output from the DAQ system, which was then fed to the EM shaker. • A resistor box which was controlled manually to apply purely resistive electrical loads in the range 0.5-500 kΩ.
The signal delivered to the shaker was bandlimited white noise with frequency bandwidth set to 0-300 Hz. In order to obtain transmissibility and voltage FRFs (frequency response functions), the accelerations at the base and top of the energy harvesting device in Figure 5, and piezoelectric voltage need to be measured. The base accelerometer (a 1 ) is connected to channel no. 2 (CH2) of the DAQ system and is set up as the reference signal (its output provides the denominator of the transmissibility and voltage FRFs). The top accelerometer (a 2 ) is connected to CH4. A resistance box is connected across the piezoelectric elements and the voltage difference across the resistance is fed to CH3. The required FRFs were generated by the DAQ system spectral analysis software from the acceleration and voltage responses channeled to the DAQ system.
. It should also be noted that, hough the DSM undamped natural frequencies are determined from the FRF, the DSM still regarded as an analytical method since the FRF terms are based on the analytical lution (Equation (17)) of the dynamic equations of the curved beam segment.

Experimental Testing
In this section, the experimental setup is presented (4.1) followed by a descriptio the estimation of the damping in the energy harvesting device (4.2). Figure 5 shows an annotated photograph of the energy harvesting device and Fig  6. shows the experimental setup. The experimental setup was comprised of the follow equipment.

Experimental Setup
• A PC-controlled data acquisition (DAQ) system (not shown in Figure 6) consistin hardware (LMS Scadas 5) and spectral analysis software (LMS Test.Lab).

•
Two accelerometers (PCB 352C22), one attached to the base of the energy harvest device and the other attached to its tip mass.

•
A signal conditioner/amplifier, to remove the noise from signals received from ac erometers and amplify them before channeling to the DAQ system.

•
An electromagnetic (EM) shaker on which the energy harvesting device w mounted.

•
A power amplifier to control the gain of the random excitation signal output fr the DAQ system, which was then fed to the EM shaker. • A resistor box which was controlled manually to apply purely resistive electr loads in the range 0.5-500 kΩ.
The signal delivered to the shaker was bandlimited white noise with frequency ba width set to 0-300 Hz. In order to obtain transmissibility and voltage FRFs (freque response functions), the accelerations at the base and top of the energy harvesting dev in Figure 5, and piezoelectric voltage need to be measured. The base accelerometer (a connected to channel no. 2 (CH2) of the DAQ system and is set up as the reference sig (its output provides the denominator of the transmissibility and voltage FRFs). The accelerometer (a ) is connected to CH4. A resistance box is connected across the piez lectric elements and the voltage difference across the resistance is fed to CH3. The requi FRFs were generated by the DAQ system spectral analysis software from the accelerat and voltage responses channeled to the DAQ system.

Estimation of Damping
To estimate the damping ratio from the experimental data, we approximate the device as a base-excited single-degree-of-freedom system, for which the transmissibility can be expressed as a function of excitation frequency , the undamped natural frequency , and viscous damping ratio :

Estimation of Damping
To estimate the damping ratio from the experimental data, we approximate the device as a base-excited single-degree-of-freedom system, for which the transmissibility can be expressed as a function of excitation frequency ω, the undamped natural frequency ω n , and viscous damping ratio ζ: From Equation (37), it is seen that the expression for T − 1 = (ω/ω n ) 2 / 1 − (ω/ω n ) 2 + 2(ω/ω n )ζj is of the same form as the accelerance (inertance) FRF of a viscously damped single-degree-of-freedom system [44]. Hence, the graph of Im{T − 1} (y-axis) vs. Re{T − 1} (x-axis) (i.e., the Nyquist plot of T − 1) will be approximately circular, just like for the accelerance FRF [45]. The Nyquist plot of T will be the same as that for T − 1 but shifted by 1 unit to the right along the horizontal axis, as shown in Figure 7. From Equation (37), it is seen that the expression for − 1 = ( ⁄ ) 1 − ( ⁄ ) + 2( ⁄ ) j ⁄ is of the same form as the accelerance (inertance) FRF of a viscously damped single-degree-of-freedom system [44]. Hence, the graph of Im − 1 (y-axis) vs. Re − 1 (x-axis) (i.e., the Nyquist plot of − 1) will be approximately circular, just like for the accelerance FRF [45]. The Nyquist plot of will be the same as that for − 1 but shifted by 1 unit to the right along the horizontal axis, as shown in Figure 7. Three approaches based on the expression in Equation (37) were used to estimate the damping ratio from the experimental transmissibility data.
The first approach considered that Equation (38a) was used to locate the resonance point ( = ), thus enabling to be determined from Equation (38b).
The second approach considered that, | ( = ) − 1| = 1 (2 ) ⁄ , the quality factor [45]. Hence, the half-power points are defined as two frequency points and (assuming ) on either side of , for which Substituting Equation (37) into Equation (39), these frequencies are determined from the following equation: Figure 7. Nyquist plot of transmissibility function for single degree of freedom system (Equation (37)).
Three approaches based on the expression in Equation (37) were used to estimate the damping ratio from the experimental transmissibility data.
The first approach considered that Equation (38a) was used to locate the resonance point (ω = ω n ), thus enabling ζ to be determined from Equation (38b).
The second approach considered that, |T(ω = ω n ) − 1| = 1/(2ζ), the quality factor [45]. Hence, the half-power points are defined as two frequency points ω 1 and ω 2 (assuming ω 2 > ω 1 ) on either side of ω n , for which Substituting Equation (37) into Equation (39), these frequencies are determined from the following equation: Considering small damping ζ < 1 2 √ 2 and neglecting terms in ζ 2 and higher order terms, the standard half-power point relation is obtained: Hence, after determining the two frequencies for which Equation (40) applies, ζ is determined from Equation (41).
The third approach for estimating the damping was to perform a circle-fit of the Nyquist plot of the experimental transmissibility data. Similar to the Nyquist plots of the receptance/mobility accelerance FRFs of viscously damped single-degree-of-freedom systems [44], the resonance point ω = ω n was identified by the location where the angular spacing of the data points was a maximum and the half power points ω = ω 1 , ω 2 subsequently determined by displacing 90 • on either side of the resonance frequency location, as shown in Figure 7.
The damping values obtained using these three methods are presented in Table 4 for different magnitudes of the external electrical load. The damping ratio estimates by the three methods are in good agreement, however the circle-fit method is expected to be the most accurate since it is less prone to errors introduced by the frequency resolution of the data and measurement noise.   With reference to the results in Table 4, for any given method, the damping ratios are minimal at the lowest resistance (approximate short circuit) and highest resistance (approximate open circuit). In these cases, the damping ratio ζ is due to the mechanical effects only since very little electrical power is dissipated due to the resistance being very low (short circuit) or the current being very low (open circuit). In between these two extremes, the estimated damping ratio is higher since it contains the additional damping effect introduced by the energy harvesting. It should also be noted that the undamped natural frequency ω n in Equation (37) will vary with the magnitude of the electrical load (since the electrical load affects the effective stiffness, apart from the damping). The transmissibility resonance condition was defined above as ω = ω n . For 0 < ζ < 0.1 (as in the short and open circuit conditions, Table 4) this frequency is very close (but not identical to) the frequency at which the transmissibility is maximum. However, given that the energy harvesting effect raises ζ to a value close to 0.1, it is important to make the distinction between the two frequencies. Hence, for the remainder of the paper the following terms are used to distinguish from the condition ω = ω n : • Transmissibility resonance peak frequency-frequency at which transmissibility FRF magnitude is maximum; • Voltage FRF resonance peak frequency-frequency at which the voltage FRF magnitude is maximum.
It is also noted that, in the analysis results presented in the following section, the value used for β (= 2ζ 1 /ω 1 , Equation (8)) uses values of ζ 1 and ω 1 obtained from experimental data for the short circuit condition (lowest resistance in Table 4). This value is used in the DSM, as well as in the FE as the stiffness-proportional damping multiplier ([C] = β[K]).

DSM vs. 2D/3D ANSYS vs. Experiment
A harmonic response analysis was performed in ANSYS to obtain the transmissibility FRF using the different types of elements. With reference to Figures 1 and 3, the boundary condition for the x-displacement amplitude of the base was set to unity and the x-displacement amplitude of the tip was therefore considered as the transmissibility. Figure 8 shows the predicted transmissibility FRFs by the different FE models, together with DSM prediction given by Equation (31). These predictions are for the case of zero electrical load (pure short circuit). The experimental result for the lowest electrical load considered (0.5 kΩ) is also included; although the experimental load is not zero, it is low enough to be considered as "short circuit" in practice. The resonance frequencies of the plane, shell, and solid element model transmissibility results in Figure 8 are very close to those already given in Table 3 since the system is relying only on the mechanical damping and the resonance frequency is weakly affected by damping if the damping ratio is less than 0.1. The result from the DSM model agrees very closely with that from the plane element model, as expected since both models describe a beam-bending behavior.
Sensors 2022, 22, x FOR PEER REVIEW 17 • Transmissibility resonance peak frequency-frequency at which transmissib FRF magnitude is maximum; • Voltage FRF resonance peak frequency-frequency at which the voltage FRF ma tude is maximum.
It is also noted that, in the analysis results presented in the following section, value used for (= 2 / , Equation (8)) uses values of and obtained from perimental data for the short circuit condition (lowest resistance in Table 4). This valu used in the DSM, as well as in the FE as the stiffness-proportional damping multip ( = ).

DSM vs. 2D/3D ANSYS vs. Experiment
A harmonic response analysis was performed in ANSYS to obtain the transmissib FRF using the different types of elements. With reference to Figure 3 and Figure 1 boundary condition for the x-displacement amplitude of the base was set to unity and x-displacement amplitude of the tip was therefore considered as the transmissibility. ure 8 shows the predicted transmissibility FRFs by the different FE models, together w DSM prediction given by Equation (31). These predictions are for the case of zero elect load (pure short circuit). The experimental result for the lowest electrical load consid (0.5 kΩ) is also included; although the experimental load is not zero, it is low enoug be considered as "short circuit" in practice. The resonance frequencies of the plane, s and solid element model transmissibility results in Figure 8 are very close to those alre given in Table 3 since the system is relying only on the mechanical damping and the onance frequency is weakly affected by damping if the damping ratio is less than 0.1. result from the DSM model agrees very closely with that from the plane element mo as expected since both models describe a beam-bending behavior. Next, the piezoelectric effect is introduced into the models by applying a signifi magnitude of electrical resistance. Figure 9 shows the predicted transmissibility and v age FRF using three different models for the case of 50 kΩ resistance, together with corresponding experimental results. The peaks of the transmissibility FRFs in Figure 9  Next, the piezoelectric effect is introduced into the models by applying a significant magnitude of electrical resistance. Figure 9 shows the predicted transmissibility and voltage FRF using three different models for the case of 50 kΩ resistance, together with the corresponding experimental results. The peaks of the transmissibility FRFs in Figure 9 are all seen to be lower than those in Figure 8, due to the additional damping introduced by the piezoelectric effect. There is again very close agreement between the DSM and plane element model predictions, both with regard to transmissibility and voltage FRF (the DSM voltage FRF being computed using Equation (35)).
(than the FE solid element model) since their resonance frequencies are closer to the experimental one. The likely reason why the beam-bending model resonance frequency is ~5% lower than that from the solid element model was already discussed above. The fact that the beam-bending model resonance is closer to the experimental resonance (than the FE solid element model) is somewhat unexpected but can be attributed to random uncertainties in the various modelling parameters. It should also be noted that the experiments were conducted with the device (Figure 1) vertically oriented, as shown in Figure 6, to minimize sliding friction in the linear bearing; this means that the static load from the weight of the tip mass (~4 N) slightly increased the curvature of the transducers-the increased curvature would have the effect of slightly lowering the fundamental resonance frequency of the device [38].
(a) (b) Figure 9. Predicted transmissibility (a) and voltage FRF (b) using three different models for the case of 50 kΩ resistance, together with the corresponding experimental results.  (Figures 10a,c, 11a,c, and 12a,c) it is evident that both theoretical results (DSM, ANSYS) show very close agreement. The ANSYS results presented are those from the plane element model, since these give the best agreement with DSM and experiment as shown in Section 5.1. Considering Figures 8 and 9, all models agree reasonably well with the experimental results, but the beam-bending ones (DSM, FE plane element) provide better agreement (than the FE solid element model) since their resonance frequencies are closer to the experimental one. The likely reason why the beam-bending model resonance frequency is~5% lower than that from the solid element model was already discussed above. The fact that the beam-bending model resonance is closer to the experimental resonance (than the FE solid element model) is somewhat unexpected but can be attributed to random uncertainties in the various modelling parameters. It should also be noted that the experiments were conducted with the device (Figure 1) vertically oriented, as shown in Figure 6, to minimize sliding friction in the linear bearing; this means that the static load from the weight of the tip mass (~4 N) slightly increased the curvature of the transducers-the increased curvature would have the effect of slightly lowering the fundamental resonance frequency of the device [38].

DSM and 2D ANSYS vs. Experiment over a Range of Electrical Loads
Figures 10a-d, 11a-d and 12a-d, respectively, show the voltage FRF magnitude, transmissibility magnitude, and the power FRF over a range of electrical loads. With reference to the predictions (Figures 10a,c, 11a,c and 12a,c) it is evident that both theoretical results (DSM, ANSYS) show very close agreement. The ANSYS results presented are those from the plane element model, since these give the best agreement with DSM and experiment as shown in Section 5.1.   quency is higher than that produced for excitation at the open circuit resonance peak frequency since the voltage resonance peak frequency is close to its short circuit condition. After exceeding this particular value of load, the voltage resonance peak frequency shifts to a value close to its open circuit condition [24]. After this switch, the output voltage magnitude is less sensitive to variation in load (as seen from Figure 10c,d, where the curves at the highest three loads are very close to each other). With regard to the transmissibility FRFs (Figure 11a-d), the agreement between simulations and experiments is satisfactory overall. The abrupt shift in resonance peak frequency beyond ~10 kΩ observed previously in the voltage FRFs is also evident in the predicted and measured transmissibility graphs. It should be noted that the values for resonance peak frequencies predicted/observed in the voltage FRFs are slightly different from those predicted/observed in the transmissibility FRFs due to the forms of the voltage and transmissibility FRFs being different [39]. It is evident from Figure 11c,d that damping is minimal (transmissibility peak magnitude maximal) in the short circuit condition (0. Focusing on the voltage FRFs (Figure 10a-d) it can be observed that by increasing the electrical load from 0.5 to 500 kΩ, there is a noticeable shift in resonance peak frequency and a significant increase in the magnitude of the voltage in both simulated and measured data. The changes produced by increasing the load from 100 to 500 kΩ are negligible, thus the 500 kΩ case can be considered as the open circuit condition. Regarding correlation between theory and experiment for the voltage FRFs, this is satisfactory, and the discrepancies in voltage FRF magnitudes reduce with increasing electrical load. As far as the voltage FRF resonance peak frequency is concerned, the measured short circuit and open circuit resonance peak frequencies are 46.7 Hz and 48.925 Hz, respectively. The corresponding predictions are 47.06 Hz and 49.08 Hz using DSM, and 47.3 Hz and 49.2 Hz using ANSYS. For loads of up to approximately 10 kΩ it can be observed that the (predicted/measured) output voltage for excitation at the short circuit resonance peak frequency is higher than that produced for excitation at the open circuit resonance peak frequency since the voltage resonance peak frequency is close to its short circuit condition. After exceeding this particular value of load, the voltage resonance peak frequency shifts to a value close to its open circuit condition [24]. After this switch, the output voltage magnitude is less sensitive to variation in load (as seen from Figure 10c,d, where the curves at the highest three loads are very close to each other).
With regard to the transmissibility FRFs (Figure 11a-d), the agreement between simulations and experiments is satisfactory overall. The abrupt shift in resonance peak frequency beyond~10 kΩ observed previously in the voltage FRFs is also evident in the predicted and measured transmissibility graphs. It should be noted that the values for resonance peak frequencies predicted/observed in the voltage FRFs are slightly different from those predicted/observed in the transmissibility FRFs due to the forms of the voltage and transmissibility FRFs being different [39]. It is evident from Figure 11c,d that damping is minimal (transmissibility peak magnitude maximal) in the short circuit condition (0.5 kΩ) and open circuit condition (500 kΩ). The damping is maximal (transmissibility peak magnitude minimal) at~10 kΩ in the predictions and between 10 and 50 kΩ in the measurement. This is in line with the estimates of the damping ratio in Table 4. It is also in line with the power FRFs (Figure 12a-d), where it is seen that the peak power dissipated is minimal in the short and open circuit conditions (0.5 kΩ and 500 kΩ, respectively), and maximal at~10 kΩ in the predictions and between 10 and 50 kΩ in the measurement. Figure 13a,b illustrates this relation more clearly between damping and energy harvesting by showing the variation of the peak transmissibility and peak power with electrical load (more experiments are required to precisely locate the experimental resistance value for minimal peak transmissibility/maximum power dissipation).
with the power FRFs (Figure 12a-d), where it is seen that the peak power dissipate minimal in the short and open circuit conditions (0.5 kΩ and 500 kΩ, respectively), maximal at ~10 kΩ in the predictions and between 10 and 50 kΩ in the measurem Figure 13a,b illustrates this relation more clearly between damping and energy harves by showing the variation of the peak transmissibility and peak power with electrical l (more experiments are required to precisely locate the experimental resistance value minimal peak transmissibility/maximum power dissipation).
The measured variations of the peak transmissibility and peak power with elect load shown in Figure 13 are consistent with the corresponding predicted variations, there are significant differences in the predicted and measured values. These are tributed to uncertainty in the various modelling parameters input into the model, uncertainty in the mechanical and electrical characteristics of THUNDER's layer ma als, manufacturing imperfections resulting in deviations from the assumed circular sh and non-uniformity of the thickness of the layers; errors in modelling of the joints link the THUNDER transducers to the rest of the energy harvester. Energy harvesting research typically focuses on the magnitude of the FRFs, [31,36,37]. The transmissibility and voltage FRFs are complex-valued and thus have b magnitude and phase information. Rather than looking at the phase variation in isolat a more thorough approach for validating the predicted effect of the phase variation plot the predicted and measured FRFs on a complex plane (imaginary part vs. real p i.e., as Nyquist plots (as done in [24] for a simple base-excited cantilever harves Nyquist plots for theoretical and experimental FRF data are shown in Figure 14. The perimental transmissibility Nyquist plots in Figure 14b were used to generate the da ing ratio data in the last column of Table 4 (after circle fitting). From Figure 14a,b it is c that the transmissibility FRF continues to follow Equation (41) as the electrical load i creased, its Nyquist plot remaining approximately circular and the orientation of the c practically unaffected by the magnitude of the load. The Nyquist plots for the predi and measured voltage FRFs (Figure 14c,d) are also approximately circular for all elect loads but their orientation is highly dependent on the magnitude of the electrical l The satisfactory agreement between theoretical and experimental Nyquist plots for b transmissibility and voltage FRF provides a further source of validation for the anal in this paper. The measured variations of the peak transmissibility and peak power with electrical load shown in Figure 13 are consistent with the corresponding predicted variations, but there are significant differences in the predicted and measured values. These are attributed to uncertainty in the various modelling parameters input into the model, e.g., uncertainty in the mechanical and electrical characteristics of THUNDER's layer materials, manufacturing imperfections resulting in deviations from the assumed circular shape, and non-uniformity of the thickness of the layers; errors in modelling of the joints that link the THUNDER transducers to the rest of the energy harvester.
Energy harvesting research typically focuses on the magnitude of the FRFs, e.g., [31,36,37]. The transmissibility and voltage FRFs are complex-valued and thus have both magnitude and phase information. Rather than looking at the phase variation in isolation, a more thorough approach for validating the predicted effect of the phase variation is to plot the predicted and measured FRFs on a complex plane (imaginary part vs. real part), i.e., as Nyquist plots (as done in [24] for a simple base-excited cantilever harvester). Nyquist plots for theoretical and experimental FRF data are shown in Figure 14. The experimental transmissibility Nyquist plots in Figure 14b were used to generate the damping ratio data in the last column of Table 4 (after circle fitting). From Figure 14a,b it is clear that the transmissibility FRF continues to follow Equation (41) as the electrical load is increased, its Nyquist plot remaining approximately circular and the orientation of the circle practically unaffected by the magnitude of the load. The Nyquist plots for the predicted and measured voltage FRFs (Figure 14c,d) are also approximately circular for all electrical loads but their orientation is highly dependent on the magnitude of the electrical load. The satisfactory agreement between theoretical and experimental Nyquist plots for both transmissibility and voltage FRF provides a further source of validation for the analysis in this paper.

Conclusions
This paper has presented an analytical model of a piezoelectric energy harvesting curved beam based on the dynamic stiffness method (DSM) and applied it to predict the measured output of an energy harvester that uses a commercial curved transducer (THUNDER TH-7R). The design of the device is novel for energy harvesting purposes in that the input excitation is applied to one end of the parallel transducers in the longitudinal direction, vibrating a tip mass attached at the other end. The DSM modelled the transducer as a curved laminated beam, with no limitations on the degree of initial curvature, and including all layers of the transducer. As a secondary novel contribution, the same device was also modelled using commercial FE software (ANSYS) through alternative element-type models. The results of DSM were consistent with those of ANSYS. Predictions by both ANSYS and DSM for the transmissibility and voltage FRFs at various electrical loads ranging from short circuit to open circuit showed a satisfactory degree of correlation with the experimental results in all aspects: magnitude, phase, Nyquist plots, and resonance frequency shift. The range of electrical load for maximal power generation and maximal damping was correctly identified. The shift in resonance frequency while changing the electrical load from short circuit to open circuit was predicted to be 4.01% by ANSYS and 4.3% by DSM, which compared well with the shift of 4.7% observed from the experimental data.
It is concluded that the presented dynamic stiffness model is sufficiently accurate for commercial curved transducers used in applications where the predominate mode of vibration is beam-like. In such situations, it can be used in preference to FE software with less computational effort due to the drastically reduced number of elements required and

Conclusions
This paper has presented an analytical model of a piezoelectric energy harvesting curved beam based on the dynamic stiffness method (DSM) and applied it to predict the measured output of an energy harvester that uses a commercial curved transducer (THUNDER TH-7R). The design of the device is novel for energy harvesting purposes in that the input excitation is applied to one end of the parallel transducers in the longitudinal direction, vibrating a tip mass attached at the other end. The DSM modelled the transducer as a curved laminated beam, with no limitations on the degree of initial curvature, and including all layers of the transducer. As a secondary novel contribution, the same device was also modelled using commercial FE software (ANSYS) through alternative elementtype models. The results of DSM were consistent with those of ANSYS. Predictions by both ANSYS and DSM for the transmissibility and voltage FRFs at various electrical loads ranging from short circuit to open circuit showed a satisfactory degree of correlation with the experimental results in all aspects: magnitude, phase, Nyquist plots, and resonance frequency shift. The range of electrical load for maximal power generation and maximal damping was correctly identified. The shift in resonance frequency while changing the electrical load from short circuit to open circuit was predicted to be 4.01% by ANSYS and 4.3% by DSM, which compared well with the shift of 4.7% observed from the experimental data.
It is concluded that the presented dynamic stiffness model is sufficiently accurate for commercial curved transducers used in applications where the predominate mode of vibration is beam-like. In such situations, it can be used in preference to FE software with less computational effort due to the drastically reduced number of elements required and the ease of application on general purpose coding software, such as MATLAB. Moreover, the matrix assembly procedure presented allows application to more complex transducers comprising segments of different radius of curvature.
Funding: This research received no external funding.

Institutional Review Board Statement: Not applicable.
Informed Consent Statement: Not applicable.

Data Availability Statement:
The data presented in this study can be made available upon request to the corresponding author.

Conflicts of Interest:
The authors declare no conflict of interest.