*Article* **Numerical Analysis and Experimental Verification of Damage Identification Metrics for Smart Beam with MFC Elements to Support Structural Health Monitoring**

**Andrzej Koszewnik 1,\*, Kacper Lesniewski <sup>1</sup> and Vikram Pakrashi <sup>2</sup>**


**Abstract:** This paper investigates damage identification metrics and their performance using a cantilever beam with a piezoelectric harvester for Structural Health Monitoring. In order to do this, the vibrations of three different beam structures are monitored in a controlled manner via two piezoelectric energy harvesters (PEH) located in two different positions. One of the beams is an undamaged structure recognized as reference structure, while the other two are beam structures with simulated damage in form of drilling holes. Subsequently, five different damage identification metrics for detecting damage localization and extent are investigated in this paper. Overall, each computational model has been designed on the basis of the modified First Order Shear Theory (FOST), considering an MFC element consisting homogenized materials in the piezoelectric fiber layer. Frequency response functions are established and five damage metrics are assessed, three of which are relevant for damage localization and the other two for damage extent. Experiments carried out on the lab stand for damage structure with control damage by using a modal hammer allowed to verify numerical results and values of particular damage metrics. In the effect, it is expected that the proposed method will be relevant for a wide range of application sectors, as well as useful for the evolving composite industry.

**Keywords:** Structural Health Monitoring; piezoelectric; energy harvesting; damage detection; macro fiber composites (MFC); damage sensitive feature; finite element method (FEM)

#### **1. Introduction**

Modern mechanical and civil structures are becoming increasingly flexible and complex with time. Beams, pipes, and cables in disparate areas of engineering (e.g., aeronautical, bridge, petroleum, renewable energy) remain critical structural elements which continue to degrade and remain susceptible to both external excitation and internal disturbances, leading to increasing risk and maintenance costs [1–5]. Critical components in these complex structures are typically designed around limit state principles for a certain level of damage tolerance and their maintenance and inspection schedules may not necessarily provide an appropriate evolution of damage due to logistic and procedural aspects, epistemic and aleatory uncertainties around such processes, and due to human effects.

In this context, damage detection or assessment of the current condition of structure can be addressed effectively through Structural Health Monitoring (SHM) by reducing the number and frequency of inspections, uncertainties, and also provide some information about the capacity of a system, which is not possible from the more popular visual inspections, thereby also increasing the value of information from such systems and eventually create more resilient structures. A wide range of SHM systems are widely studied by both academia and the industry, including vibration, optical, thermal, or impedance-based

**Citation:** Koszewnik, A.; Lesniewski, K.; Pakrashi, V. Numerical Analysis and Experimental Verification of Damage Identification Metrics for Smart Beam with MFC Elements to Support Structural Health Monitoring. *Sensors* **2021**, *21*, 6796. https://doi.org/10.3390/s21206796

Academic Editor: Zdenek Hadas

Received: 19 August 2021 Accepted: 8 October 2021 Published: 13 October 2021

**Publisher's Note:** MDPI stays neutral with regard to jurisdictional claims in published maps and institutional affiliations.

**Copyright:** © 2021 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (https:// creativecommons.org/licenses/by/ 4.0/).

methods [6–8]. In particular, the vibration method is intensively developed by many researchers due to the simplicity of implementation of chosen sensors on real structures and fast detection of damage in the structures with the use of real-time monitoring systems. The instance of these considerations are papers [9,10], where the behavior of the intact and damaged mechanical structures has been assessed based upon natural frequencies, mode shapes, frequency response functions, and also power spectral density or based on spatial wavelet analysis [11].

Many investigations in this field have shown that one-layer piezoelectric patches or piezo harvesters with lightweight fiber materials allow identifying changes in different kinds of vibrating structures [12–16]. The relevant properties of these element, especially their fragility or extreme value effects allow to use them in many applications to monitor structure or detect damage [17–21]. A core aspect behind energy harvesting not being a part of SHM is the relative lack of sector specific examples and benchmarks, where the damage detection markers are compared in terms of their performance. Such markers can relate to damage existence, location, extent, and a combination thereof and it is important to discover and distinguish which markers are relevant for which purposes. In recent times, there has been some effort in creating benchmarks yielding sector specific challenges and energy harvesting sensors, along with markers for damage detection. A recent work [22] uses energy harvesting (EH) systems to assess the leak localization in water pipes. Here, the authors investigated several pipes with different widths of leak to propose and calibrate a leak index based on the monitoring voltage from a piezoelectric energy harvester (PEH) and the power spectrum of the output signal generated from particular polyvinylidene fluoride (PVDF) piezoelectric transducers. Subsequently, the use of Pb-free biomolecular piezoelectrics was also used to enhance SHM of water pipes [23]. Similar examples are available for bridge monitoring under operational conditions via Pb-zirconate titanate (PZT) patches [24,25]. Similar investigations have been also performed with piezoelectric PVDF sensors for wireless monitoring of tension conditions in a cable stayed bridge [26]. Similar applications have also been considered for the aerospace sector in terms of component monitoring in airplanes [27], including those harvesting energy from fuselage vibration with one-layer piezo transducers [28].

While the potential of an energy harvesting system for supporting SHM is established in principle, the implementation of it in specific engineering sectors is still fraught with several questions around interpretation and performance of responses, their analyses, and metrics developed for monitoring features of interest. An initial computational model is helpful in this regard to assess sensor location, potential impact of nonlinearities of the energy harvesting element, and the eventual translation of such information into designing appropriate SHM systems. The advantage of an energy harvesting-based monitoring often lies in the low-power aspect of it and the ease of use, which can lead to the possibility of extensive instrumentation. Model updating and digital twinning are also becoming more common in such industries and, consequently, a work like this will also provide a connectivity of harvesters integrated to such updating processes, for operational structures and those which are evolving through varied technological readiness levels.

There exists papers which address elements of the abovementioned challenges. The influence of non-linear geometric responses of a piezoelectric composite plate considering von Karman large strain theories into the classic plate solution has been investigated by using a 3D element model [29], with results indicating that the problem cannot be omitted especially when correct prediction of the stress-strain over the PEH is analyzed. A piezoelectric element modeled as a shell element acting under d31 effect of the crystal [30] noted that the effect of non-linearity is small and can be neglected, especially when commercial piezoelectric patches are used. Other scientific works in this field focused on introducing piezoelectric coupling to the shell element [31,32] and the results indicate the influence of non-linearity for piezoelectric laminated shell is significant and should be further analyzed. As a result, this led to developing investigations and modeling the piezo structure as a higher order layer-wise plate finite element considering piezoelectric coupling [33]. In

summary, despite many scientific contributions related to formulations of plates and shells for piezoelectric laminated elements, there is a gap in verification of numerical results considering the shell finite elements of the piezoelectric element. This paper addresses this knowledge gap by carrying out numerical analysis and subsequently validating them against experimental results.

Studies similar to what has been presented in this paper are also of particular relevance for new sensors that are being developed from environmental perspectives to avoid Pbbased systems [23] or multi-functional materials [34]. As their material properties and uncertainties become lower, the possibility of their use in many sectors get higher and the current work can make them better adapted to the composites sector where energy harvesting based SHM still requires significant research. While the presence of damage is often easier to detect with a sensor, the localization of such damages often provides further challenges, especially when the damages are relatively closely spaced. The impact of closely spaced damages have been studied before [11,20], and while there is a natural understanding that indicators of damage at some point would present a coalesced effect of two damages after a certain proximity, there is a paucity of literature in investigating quantitatively what the effects are for patch based systems. Furthermore, the use of composite materials for this topic is important as a benchmark due to their extensive use in new sectors, especially in marine environments [21]. Over time, their degradation, especially in the durability aspects for saline and harsh marine environments [35], will be of particular relevance around this topic. There is thus a need for a detailed numerical and experimental investigation of a relatively generic example which can be used in the future for similar studies, but also as an evidence base for current performance of patchbased energy harvesting and SHM, future adaption to new sensors, and to new structural systems and environments. This will also lend eventually to estimates on their lifetime risk levels [36] and comparative performances [37] around such levels, especially when the exposure is of a stochastic nature.

Therefore, this paper is focused on damage identification in thin structures by using piezo fiber composites. In contrast to [20,21,23], the piezo harvester with three dimensional material in the piezoelectric fiber layer is considered to enhance the harvesting effect, especially in resonance regions with the host structure. To establish this, three different structures have been investigated, where the first one is intact, while the other two are damaged. The damage is introduced by drilling holes in the area of the beam closer to the end of the piezo-composite located closer the fixed end of the beam. The numerical results obtained from the finite element (FE) models of both sensors and also the damage indexes were determined next based on the frequency response functions, which subsequently allows for damage localization assessment. An experimental test on a real structure was finally carried out to verify the numerical results.

The paper is organized as follows: Section 2 describes the methodology of modeling a piezo-harvester based on the First Order Shear Theory (FOST) and also presents the procedure of determining the shell finite element of the laminated structure of the piezocomposite. Section 3 presents the computational models of all structures (intact and damage) with a homogenous model of macro-fiber composite (MFC), which is also a core novelty of this manuscript. In Section 4, experimental investigations were carried out for all structures with the use of an impact modal hammer. The measurement signals from both MFC elements allow calculation of real values of damage matrices and for comparison with numerical results. Additionally, the approach allows assessing the damage location. Section 5 concludes the main findings of this work.

#### **2. Mathematical Model of the Smart Structure with a Laminated Model of the MFC Element**

The finite element model designed on the basis of the First Order Shear Theory (FOST) is mainly used in many applications consisting of thin structures like beams or plates to monitor their state via piezoelectric patches like macro-fiber composites [38]. The chosen finite element (FE) model is described by the laminated shell quadratic finite element with eight nodes: five mechanical degrees of freedom (three translations, two rotations) and three electrical DOFs related to the number of the piezoelectric-active layers of MFCs. In addition, this FE model can act as an equivalent single layer model to describe the mechanical behavior and as a layer-wise model to describe the electrical behavior. As a result, all the aforementioned performances of the proposed FE element caused that this element can simulate not only layered materials (e.g., laminated structural composites) with the piezoelectric layers embedded in the structure, but also non-layered materials (e.g., aluminum) with piezoelectric sensors attached to this structure. Another important aspect to note is that only displacements, forces applied to the structure (impulse inputs), and electrical potentials can be enough to obtain the frequency response functions of this structure and calculate the damage metrics properly. This makes the proposed FE model more useful to increase the computation efficiency, especially for real time SHM, without compromising damage identification accuracy.

The shell finite element, capable of simulating a plane, can be implemented in Ansys and solved using the curvature formulation, which must be at least quadratic to describe single or double curved shells with adequate accuracy [39]. Under such circumstances, the problem of interpolation of the curvature field should be considered in the same way as in the case of the issue of variables. In order to do this, it can be assumed that direction 1 is taken to be aligned to the fibers of a given lamina, direction 3 is aligned to the normal laminate direction, while direction 2 is obtained based on the right-hand rule. Indexes of all variables described in this manuscript are then given in the range of *i* = 1..3, *j* = 1..3, *k* = 1..3, *l* = 1..3.

#### *Constitutive Equations and Electrical Assumptions of the Piezo-Electric Lamina and Its Finite Element Model*

The piezoelectric MFC, as shown in Figure 1, is a five-layer smart composition with a single active layer of the thickness of *tMFCa*, two electrode layers of the thickness of *tMFCe,* and two Kapton layers of the thickness of *tMFCk*. As a result, this smart MFC element, attached to the host structures, is described in the form of constitutive equations given by Equation (1).

$$\begin{aligned} \sigma\_{ij} &= \mathbb{C}^{E}\_{ijkl}\varepsilon\_{kl} - \mathfrak{e}\_{kij}\mathrm{E}\_{k\prime} \\ D\_{i} &= \mathfrak{e}\_{ijk}\varepsilon\_{kl} + d^{\varepsilon}\_{ik}\mathrm{E}\_{k} \end{aligned} \tag{1}$$

where:

*σij*—the laminate stresses;

*Dk*—the electrical displacement;

*Ek*—the electrical field;

*CE ijkl*—the short-circuit elastic properties of the piezo-laminate;

*εkl*—the laminate strains;

*eijk*—the piezoelectric coupling coefficients;

*dε ik*—the dielectric properties of the piezo-laminate.

**Figure 1.** The structure of the macro-fiber composite.

From the measurement point of view, a proper polarization direction of the active element attached to the structure is required. In many applications, especially in those related to SHM, this problem has been solved by choosing the transverse direction of the active element polarization. As a result, Equation (1) can be simplified to the form given by Equation (2) when the Voigt notation is used, and also the normal transversal stress can be considered by omitting (*σ*<sup>33</sup> ≈ 0).

$$\begin{Bmatrix} \sigma\_{11} \\ \sigma\_{22} \\ \sigma\_{12} \\ \tau\_{23} \\ \tau\_{13} \\ D\_3 \end{Bmatrix} = \begin{bmatrix} Q\_{11} & Q\_{12} & 0 & 0 & 0 & -\varepsilon'\_{31} \\ Q\_{12} & Q\_{22} & 0 & 0 & 0 & -\varepsilon'\_{31} \\ 0 & 0 & Q\_{66} & 0 & 0 & 0 \\ 0 & 0 & 0 & Q\_{44} & 0 & 0 \\ 0 & 0 & 0 & 0 & Q\_{55} & 0 \\ \varepsilon'\_{31} & \varepsilon'\_{31} & 0 & 0 & 0 & d'\_{33} \end{bmatrix} \begin{Bmatrix} \varepsilon\_{11} \\ \varepsilon\_{22} \\ \varepsilon\_{12} \\ \gamma\_{23} \\ \gamma\_{31} \\ E\_3 \end{Bmatrix}, \tag{2}$$

where the material matrix components are

$$\begin{split} Q\_{11} &= \mathbb{C}\_{11}^{\mathrm{E}} - \frac{\mathbb{C}\_{13}^{\mathrm{E}2}}{\mathbb{C}\_{33}^{\mathrm{E}}}, Q\_{12} = \mathbb{C}\_{12}^{\mathrm{E}} - \frac{\mathbb{C}\_{13}^{\mathrm{E}2}}{\mathbb{C}\_{33}^{\mathrm{E}}}, Q\_{22} = \mathbb{C}\_{22}^{\mathrm{E}} - \frac{\mathbb{C}\_{13}^{\mathrm{E}2}}{\mathbb{C}\_{33}^{\mathrm{E}}}, Q\_{66} = 2\left(\mathbb{C}\_{11}^{\mathrm{E}} - \mathbb{C}\_{22}^{\mathrm{E}}\right), \\ Q\_{44} &= \mathbb{C}\_{44}^{\mathrm{E}} + \frac{\mathbb{c}\_{15}^{2}}{\mathbb{d}\_{11}}, Q\_{55} = \mathbb{C}\_{55}^{\mathrm{E}} + \frac{\mathbb{c}\_{15}^{2}}{\mathbb{d}\_{11}}, \; \mathfrak{e}\_{31} = \mathfrak{e}\_{31} - \frac{\mathbb{C}\_{13}^{\mathrm{E}}\mathfrak{e}\_{33}}{\mathbb{C}\_{33}^{\mathrm{E}}}, \; \mathfrak{d}\_{33}^{\mathrm{d}} = \mathfrak{d}\_{33}^{\mathrm{E}} - \frac{\mathfrak{e}\_{31}^{2}}{\mathbb{d}\_{13}^{\mathrm{E}}}. \end{split}$$

Application of smart elements in the form of macro fiber composites in mechanical structures requires also considering their mechanical and electrical properties. In the case of the mechanical behavior of the laminated structures, the Equivalent Single-Layer (ESL) approach has been used. Then, according to Figure 1, actuating and generalized forces and moments, like: shearing forces (Q), bending moments (M), normal moments (N), and torsion moments (T) in relation to the thickness of the piezo, can be written in the local coordinate system in the following forms:

$$\left\{ \begin{array}{c} N\_{\rm x} \\ N\_{y} \\ N\_{\rm xy} \end{array} \right\} = \int\_{-t\_{\rm MCx}/2}^{t\_{\rm MCx}/2} \left\{ \begin{array}{c} \sigma\_{\rm xx} \\ \sigma\_{\rm yy} \\ \sigma\_{\rm xy} \end{array} \right\} dz + \sum\_{\rm op=1}^{\rm nl\_{p}} \int\_{-t\_{\rm MCx}/2}^{t\_{\rm MCx}/2} \left\{ \begin{array}{c} \sigma\_{\rm xx} \\ \sigma\_{\rm yy} \\ \sigma\_{\rm xy} \end{array} \right\} dz + \sum\_{\rm ok=1}^{\rm nl\_{k}} \int\_{-t\_{\rm MCx}/2}^{t\_{\rm MCx}/2} \left\{ \begin{array}{c} \sigma\_{\rm xx} \\ \sigma\_{\rm yy} \\ \sigma\_{\rm xy} \end{array} \right\} dz,\tag{3}$$

$$\left\{ \begin{array}{c} M\_{\mathcal{X}} \\ M\_{\mathcal{Y}} \\ M\_{\mathcal{X}\mathcal{Y}} \end{array} \right\} = \int\_{-t\_{\mathrm{MCE}/2}/2}^{t\_{\mathrm{MCE}/2}/2} z \left\{ \begin{array}{c} \sigma\_{\mathrm{xx}} \\ \sigma\_{\mathrm{yy}} \\ \sigma\_{\mathrm{xy}} \end{array} \right\} dz + \sum\_{o=1}^{nl\_{\mathrm{p}}} \int\_{-t\_{\mathrm{MCE}/2}/2}^{t\_{\mathrm{MCE}/2}/2} z \left\{ \begin{array}{c} \sigma\_{\mathrm{xx}} \\ \sigma\_{\mathrm{yy}} \\ \sigma\_{\mathrm{xy}} \end{array} \right\} dz + \sum\_{o=1}^{nl\_{\mathrm{k}}} \int\_{-t\_{\mathrm{MCE}/2}/2}^{t\_{\mathrm{MCE}/2}/2} z \left\{ \begin{array}{c} \sigma\_{\mathrm{xx}} \\ \sigma\_{\mathrm{yy}} \\ \sigma\_{\mathrm{xy}} \end{array} \right\} dz,\tag{4}$$

 *Qx Qy* = *tMFCa*/2 −*tMFCa*/2 5 6 <sup>1</sup> <sup>−</sup> <sup>4</sup>*s*<sup>3</sup> 2 *tMFCa*<sup>2</sup> *τyz τxz dz* + *nlp* ∑ *op*=1 *tMFCp*/2 −*tMFCp*/2 5 6 <sup>1</sup> <sup>−</sup> <sup>4</sup>*s*<sup>3</sup> 2 *ttotal*\_*MFCp*<sup>2</sup> *τyz τxz dz* <sup>+</sup> *nlk* ∑ *ok*=1 *tMFCk*/2 −*tMFCk*/2 5 6 <sup>1</sup> <sup>−</sup> <sup>4</sup>*s*<sup>3</sup> 2 *ttotal*\_*MFCk*<sup>2</sup> *τyz τxz dz*, (5)

 *Trx Try* = *tMFCa*/2 −*tMFCa*/2 5 6 *z* <sup>1</sup> <sup>−</sup> <sup>4</sup>*s*<sup>3</sup> 2 *tMFCa*<sup>2</sup> *τyz τxz dz* + *nlp* ∑ *op*=1 *tMFCp*/2 −*tMFCp*/2 5 6 *z* <sup>1</sup> <sup>−</sup> <sup>4</sup>*s*<sup>3</sup> 2 *ttotal*\_*MFCp*<sup>2</sup> *τyz τxz dz* <sup>+</sup> *nlk* ∑ *ok*=1 *tMFCk*/2 −*tMFCk*/2 5 6 *z* <sup>1</sup> <sup>−</sup> <sup>4</sup>*s*<sup>3</sup> 2 *ttotal*\_*MFCk*<sup>2</sup> *τyz τxz dz*, (6)

where:

*nlp*—the amount of the electrode layers of the MFC element, *op = 1..nlp*; *nlk*—the amount of the Kapton layers of the MFC element, *ok = 1..nlk*; *nl*—the total amount of layers of the MFC element: *nl = nlp + nlk +* 1; *tMFCa*—the thickness of the active layer of the MFC element; *tMFCp*—the thickness of the passive layer of the MFC element; *tMFCk*—the thickness of the Kaption layer of the MFC element; *ttotal\_MFCp*—the total thickness of the passive layers of the MFC element; *ttotal\_MFCk*—the total thickness of the Kaption layers of the MFC element.

In the case of determining the electrical behavior of this laminate piezo structure, the layer-wise approach has been used. Then, according to this method, electrical displacement of this piezo-composite for the active piezoelectric layer is expressed in the following form:

$$D\mathfrak{z} = \int\_{-t\_{M\text{F}a}/2}^{-t\_{M\text{F}a}/2} \{E\mathfrak{z}\} dz\_{\prime} \tag{7}$$

where:

*D*3—the electrical displacement of the active layer of the MFC element

*E*3—the electrical field of the MFC element with the vertical polarization.

Determining the representative mechanical behavior of the smart structure with the attached laminate to its surface is also required and is relevant for SHM. To do this, the degenerated shell theory with an implicit curvature [38] is used. Displacements, strains, and the electrical field can be written then as a function of the nodal degree of freedom of the finite element in the following form (Equations (8)–(10)), respectively:

$$
\mu\_i = \mu\_i^n \phi\_{n\prime} \ \theta\_k = \theta\_k^n \phi\_{n\prime} \ \varphi\_p = \varphi\_p^n \phi\_{n\prime} \tag{8}
$$

$$\begin{Bmatrix} \varepsilon\_{11} \\ \varepsilon\_{22} \\ \varepsilon\_{12} \\ \gamma\_{23} \\ \gamma\_{13} \end{Bmatrix} = \left( \begin{array}{c} B\_m + B\_{b0} + s\_3 B\_{b1} \\ B\_s + B\_{t0} + s\_3 B\_{t1} \end{array} \right) \widetilde{u}\_{\prime} \tag{9}$$

$$E\_{\mathfrak{B}} = B\_{\mathfrak{G}} \widetilde{\mathfrak{p}}\_{\prime} \tag{10}$$

where: *φ<sup>n</sup>* is the shape function for n-th node of the finite element, *Bm*, *Bb*0, *Bb*1, *Bs*, *Bt*0, *Bt*<sup>1</sup> are the curvature-displacement components calculated versus in-plane membrane strains (*m*), *b*<sup>0</sup> is the uniform term of in-plane bending strains, *b*<sup>1</sup> is the linear term of in-plane bending strains, *s* is the out-of-plane shearing distortions, *t*<sup>0</sup> is a uniform term of out-of-plane torsions, and *t*<sup>1</sup> is a linear term of out-of-plane torsions, respectively.

Taking in Equations (8)–(10) to account, the strains, displacements, and electrical potentials of the laminated elements can be expressed in terms of the nodal variables. Subsequently, taking a solution of the elemental equilibrium equation adopted from [38] into account, equations for the piezoelectric problem of laminate structure where *wi* and *wj* are Gauss' Quadrature weights, can be expressed as

$$\begin{aligned} \mathbf{M}\ddot{u} + \mathbf{C}\dot{u} + \mathbf{K}\_{uu}u + \mathbf{K}\_{u\rho}\boldsymbol{\varrho} &= \mathbf{F} \\ \mathbf{K}\_{\boldsymbol{\varrho}u}u + \mathbf{K}\_{\boldsymbol{\varrho}\boldsymbol{\varrho}}\boldsymbol{\varrho} &= \mathbf{Q} \end{aligned} \tag{11}$$

where:

$$\mathbf{M} = \sum\_{i,j=1}^{j=3} \det(J^{-1}) w\_i w\_j \rho \left[ \frac{h}{2} H\_0^T H\_0 + \frac{h^2}{4} \left( H\_1^T H\_0 + H\_0^T H\_1 \right) + \frac{h^3}{12} H\_1^T H\_1 \right],\tag{12}$$

$$\mathbf{K}\_{\rm ulr} = \sum\_{i,j=1}^{j=3} \det(f^{-1}) w\_i w\_j \left( \begin{bmatrix} B\_u^{\rm ulr} \end{bmatrix}^T \begin{bmatrix} A & B \\ B & D \end{bmatrix} \begin{bmatrix} B\_u^{\rm ulr} \end{bmatrix} \right) + \sum\_{i,j=1}^{j=2} \det(f^{-1}) w\_i w\_j \left( \begin{bmatrix} B\_u^I \end{bmatrix}^T \begin{bmatrix} G & G\_h \\ G\_h & H \end{bmatrix} \begin{bmatrix} B\_u^I \end{bmatrix} \right), \tag{13}$$

$$\mathbf{K}\_{\varphi u} = \sum\_{i,j=1}^{j=3} \det(f^{-1}) w\_i w\_j \left( \begin{bmatrix} B\_{ut} \end{bmatrix}^T \begin{bmatrix} \overline{\mathfrak{E}}\_1 \\ \overline{\mathfrak{E}}\_2 \end{bmatrix} \begin{bmatrix} B\_{ut} \end{bmatrix} \right) \tag{14}$$

$$\mathbf{K}\_{\varphi\varphi} = \sum\_{i,j=1}^{j=3} \det(J^{-1}) w\_i w\_j B\_{u\varphi} \, ^T \overline{d} B\_{u\varphi} \, \tag{15}$$

$$\mathbf{C} = \alpha \mathbf{M} + \beta \mathbf{K}\_{\text{uu}\_{\prime}} \tag{16}$$

$$B\_{\rm u}^{m} = \begin{bmatrix} B\_{l0} + B\_{l0} \\ B\_{l1} \end{bmatrix}, \\ B\_{\rm u}^{t} = \begin{bmatrix} B\_{l0} + B\_{l0} \\ B\_{l1} \end{bmatrix}, \\ B\_{\rm uq} = \begin{bmatrix} B\_{l0} \\ B\_{l1} \end{bmatrix}, \\ H\_{0\_{iS(n+1)+j}} = \delta\_{ij} \phi\_{n}, \\ H\_{1\_{iS(n+1)+j}} = \delta\_{(i,j+3)} \frac{h}{2} \phi\_{n} H\_{ij}^{n}. \tag{17}$$

From a numerical point of view, the electrical and mechanical degree of freedom, as well as the generalized mass, damping and stiffness matrices need to be transformed to modal coordinates in such a way that the nodal variables for a given element can be obtained in a single vector ordered by node numbering. In order to do this, the transforming matrix R should be used. Then, a general mechanical-electrical system in modal coordinates can be express in the following form, respectively:

$$\hat{\mathbf{u}} = \left\{ \begin{array}{ccccccccc} \boldsymbol{u}\_1^1 & \boldsymbol{u}\_2^1 & \boldsymbol{u}\_3^1 & \boldsymbol{\theta}\_x^1 & \boldsymbol{\theta}\_y^1 & \boldsymbol{\rho}\_{p(1)}^1 & \dots & \boldsymbol{\rho}\_{p(last)}^1 & \dots \\ \boldsymbol{u}\_1^n & \boldsymbol{u}\_2^n & \boldsymbol{u}\_3^n & \boldsymbol{\theta}\_x^n & \boldsymbol{\theta}\_y^n & \boldsymbol{\rho}\_{p(1)}^n & \dots & \dots & \boldsymbol{\rho}\_{p(last)}^n \end{array} \right\}^T,\tag{18}$$

$$
\hat{\mathbf{u}} = \mathbf{R} \left\{ \begin{array}{c} \tilde{\boldsymbol{u}} \\ \tilde{\boldsymbol{\varphi}} \end{array} \right\} , \tag{19}
$$

$$
\hat{\mathbf{K}} = \mathbf{R} \begin{bmatrix}
\ K\_{\mu\mu} & K\_{\mu\rho} \\
\ K\_{\rho\mu} & K\_{\rho\rho}
\end{bmatrix} \mathbf{R}^T \tag{20}
$$

$$
\hat{\mathbf{M}} = \mathbf{R} \begin{bmatrix} \ \mathbf{M} & \ 0 \\ 0 & 0 \end{bmatrix} \mathbf{R}^T,\tag{21}
$$

$$\mathbf{\hat{C}} = \mathbf{R} \begin{bmatrix} \mathbf{C} & 0 \\ 0 & 0 \end{bmatrix} \mathbf{R}^T,\tag{22}$$

$$\mathbf{f} = \mathbf{R} \left\{ \begin{array}{c} F \\ Q \end{array} \right\} , \tag{23}$$

$$
\mathbf{\hat{M}\ddot{u}} + \mathbf{\hat{C}\dot{u}} + \mathbf{\hat{K}\dot{u}} = \mathbf{\hat{F}} \tag{24}
$$

#### **3. Numerical Analysis of a Smart Beam as a Laminated Structure**

The numerical calculations of the cantilever beam of a length of 380 mm, width of 31 mm, and a thickness of 1.8 mm, with two piezo-stripe elements (MFC 8528 P2), consisting of a three-dimensional homogenized material in the active layer, is described in this section. The parameters of a homogenized material for the MFC taken from [40] are collected in Table 1. To assess the values of damage identification metrics, the considered structures (intact and damage structures with first one, two and three drilled holes) were modelled in the Ansys software with the assumption that the first MFC element is located in the distance of 40 mm from the fixed end of the beam, while the second one is at a distance

of 15 mm from the free end of the beam. Taking into account the structure of the cantilever beam with MFC attached to its top surface, the host beam structure is modelled by using an 8-node coupled-brick element Solid186. For the MFC element, the electrode and Kapton layers are modelled similarly with the use of a Solid186 element, while the piezoelectric fiber layer with a homogenized material is modelled by using a Solid226-node coupled brick element. In addition, it has been assumed that modelling of the adhesive layer can be omitted due to its thickness which is less than 15 μm. Finally, this leads to determining the computational model of the considered structure shown in Figure 2 consists of 120 of finite elements Solid226 of the length of 10 mm for the cantilever beam, 102 elements of type Solid226 for each passive layer (electrode and Kapton) and 102 elements of type Solid226 for the active layer.

**Figure 2.** The numerical model of the smart intact structure with both piezo composites attached to the host structure.


**Table 1.** Material properties of homogenized MFC layer of MFC8528 P2.

In the first step of numerical calculations, the eigenvalue problem of such a modelled structure is solved by using the Ansys software. For this purpose, the behavior of three different structures (one intact and two damage structures with one hole and two holes, respectively, of the diameter of 8 mm) are compared in the selected frequency range 1–400 Hz. The example of the damage structure with two holes is shown in Figure 3. The obtained eigenvalues for each structure are listed in Table 2.


**Table 2.** Values of natural frequencies of the intact structure and damage structures.

**Figure 3.** The numerical model of the damage structure with two holes drilling in region close to the end of the piezo-patch MFC1.

Taking into account the obtained result shown in Table 2, it can be noticed that the values of the first four lowest natural frequencies of the damage structures for the increasing number of holes in the structure decreases slightly. The obtained effect is insignificant from the point of view of considering an SHM system since the decrease is of about 0.5% versus the values of eigenvalues of the intact structure. This is thus can be omitted in further analysis.

Next, a harmonic analysis for each structure is performed. To do this, the computational models of intact and damage structures with shell models of the MFC elements were excited by an impulse load of 48N magnitude, while the vertical displacement was taken from specific nodes of brick models of piezo composites MFC1 and MFC2, respectively. Taking into account the results presented in Figure 4, it can be observed that the accuracy between the frequency responses of intact and damage structures is high, especially in resonance regions, where the vibration amplitudes of the damage structure are lower than that of the intact structure (see Figure 5). This leads to difficulties related to proper identification and interpretation of damage in the structure and the calculation of damage identification metrics as damage indicators. This issue is taken up in detail in the next section.

Further analysis of the harmonic responses performed for the FEM models of the intact and damage structures indicates that the distance between the location of the measurement and excitation is the cause of obtaining two different kinds of systems from a control strategy view point, linked to collocated for MFC1 sensor and non-collocated for MFC2 sensor, respectively. An example of the frequency response for the collocated system is the upper diagram shown in Figure 4, where resonance and anti-resonance frequencies are occurring. An example of the frequency response of the non-collocated system is presented in the lower diagram in Figure 4 where the increase of the distance between the point of measurement related to the location of the piezo MFC2 and the point of impact lead to omitting the first antiresonance frequency. Consequently, a proper location of the piezo-sensor on the structure should be also considered for SHM deployment.

**Figure 4.** The frequency response function for the intact and damaged structures calculated in the selected frequency range of 1–400 Hz in terms of (**a**) the piezo MFC1 location, (**b**) the piezo MFC2 location.

**Figure 5.** The frequency response function for the intact and damaged structures calculated in the selected frequency range in terms of (**a**) the piezo MFC1 location, (**b**) the piezo MFC2 location.

#### **4. Experimental Verification**

The process of frequency response function (FRF) verification of the intact and damage structures are further used to compute damage identification metrics. Two piezo-patch sensors, MFC 8528 P2, were attached to the host aluminium beam structure with the help of an adhesive UHU Plus (Figure 6). The first piezo MFC1 is located 40 mm from the fixed end, while the second one in the distance of 15 mm from the free end of the beam. The laboratory stand is retrofitted into the modal hammer developed by Bruel and Kjaer used to excite the considered structures to vibrations. The data acquisition module PXI 4496 is used to measure and record vibrations from the beam.

Taking numerical investigations of Section 3 into account, experiments on a real structure were carried out. An impulse load with a magnitude of 48 N is applied to the structure to a chosen point located 10 mm before the piezo MFC1 and 30 mm from the fixed end of the beam while structural vibrations are measured with both laminate composites, MFC1 and MFC2, respectively.

**Figure 6.** The view of a real laboratory stand during a lab test (**a**) the smart beam with both piezocomposites MFC 8528 P2, (**b**) the data acquisition module PXI 4496a.

First, the intact structure is investigated by applying the aforementioned impulse load to the beam at 4.2 s of measurement. Taking into account the recorded voltage from both piezo elements shown in Figure 7, it can be seen that the transient response measured with MFC2 is longer than that measured with MFC1. This is due to the fact that the piezo MFC1 is located closer to the fixed end of the beam where the damping is higher. Further analysis of the intact structure behavior requires transformation of input/output signals to the frequency domain and determining two frequency responses functions (FRFs) at MFC1 and MFC2 locations, respectively, as presented in Figure 8.

**Figure 7.** The excitation signal and measurement signals measured by MFC elements during analysis of the intact structure.

**Figure 8.** The comparison of the amplitude plots of the intact structure measured by using (**a**) piezo-patch MFC1, (**b**) piezo-patch MFC2.

In Figure 8, collocated piezo MFC1 and non-collocated piezo MFC2 aspects are observed and this allows to verify FRFs determined from computational models. The frequency response for piezo MFC1 has interchangeable character of the resonance and anti-resonance frequencies, while in the case of MFC2, the generated frequency response is without the first anti-resonance frequency. Thus, the distance between the sensor and the actuator is crucial to describe the behaviour of the structure.

Further analysis of the recorded frequency responses from both piezo-sensors showed also a high convergence between them, especially in the resonance areas, where the amplitude of vibrations from tests is close to the amplitude from the numerical model. In other areas it can observe mismatch between both FRFs that are due to a heterogenous adhesion between the piezo elements and the host structure and nonlinearity of the MFC, especially in the strain-displacement relation. In the effect, the real frequency response generated on the basis of the noisy measurement signal contains additional slight amplitude peaks especially in higher frequencies. However, despite these discrepancies, from monitoring point of view, it can be still conclude that those responses properly verify the numerical responses.

In the next step, an experimental test was carried out for the damage structure with one hole drilled in the distance of 25 mm from the end of the MFC1 piezo-patch sensor and 150 mm from the place of impact. Similarly to the previous case, this structure was excited to vibration again by using a modal hammer and applying the impulse load with the magnitude of 48 N at the same point. As a result, two measurement signals from both piezo-composites were measured by a PXI module that next allow to generate two separate frequency response functions showed in Figure 9.

Observing diagrams in Figure 9, it can be noticed again that the experimental frequency response functions are close to the frequency responses obtained based on the numerical model, especially in the resonance areas of the first four natural frequencies. In this case, it can be seen that the natural frequencies of the structure measured with the help of both piezo-composites, MFC1 and MFC2, are identical with those calculated from the numerical model, while their amplitudes, especially for those measured by using the MFC1, are less convergence. This behaviour results from higher damping in the real structure than it was assumed for the numerical model. Moreover, as it was just mentioned, it is caused by heterogenous adhesion between the sensor and the host structure, nonlinearity of the piezo sensors, as well as the noisy measurement signal that is used to generate real frequency response. The mismatch is also representative of typical tests.

**Figure 9.** The comparison of the amplitude diagram of the damage structure with one hole measured by using (**a**) piezopatch MFC1, (**b**) piezo-patch MFC2.

In the same way, the structure with two holes located very close to the MFC1 sensor was investigated. In this case, the second hole was located 12 mm from the first one and 37 mm from the end of the piezo MFC1. Taking into account Figure 10, it can be observed that the amplitudes of the structure vibration on the generated FRF from the piezo MFC1 are close to the amplitudes vibrations calculated based on the numerical model. Another behaviour that can be seen in the case of the FRF analysis from the piezo MFC2 is the high convergence only in resonant areas. The main reasons of this mismatch can be attributed to a heterogenous adhesion between the bottom surface of the piezo MFC2 and the top surface of the aluminium beam, nonlinearity of the piezo-composite and noisy measurement signal. Again, despite some discrepancies between them located outside the resonance areas, the FRF generated from the lab stand can be assumed as correct.

**Figure 10.** The comparison of the amplitude plots of the damage structure with two drilled holes measured by using (**a**) piezo-patch MFC1, (**b**) piezo-patch MFC2.

The last step of the experimental test was collecting all the generated FRF from both piezo sensors—MFC1 and MFC2—to perform their analysis and assess the real value of the damage identification metrics. Taking into account Figure 11, it can be observed that the increasing number of holes in the damage structure and the decreasing stiffness of the structure in chosen areas of the beam do not lead to a change in the natural frequencies but affects only the amplitude of the structure in the resonance areas. It this way, the conclusion from the analysis of the computational model has been verified. Further analysis of these diagrams indicates also that the decrease of the beam stiffness resulting from drilling the holes in the areas located very close the MFC1 sensor leads also to the amplitude increase of vibrations measured by the MFC1 sensor but only for the first lowest natural frequencies. In the case of the piezo MFC2, it can be observed that drilling one hole leads firstly to the increase of the vibration amplitude while drilling another hole leads to its decrease. A similar effect is also presented in Figure 12 where the power spectrum of measurement signals from both piezo-sensors is analyzed. Finally, taking into account the generated FRF's and the power spectrum of signals from the piezo composites, it can be concluded that those diagrams are insufficient to identify the damage in the structure properly. For this reason, the damage identification metrics should be calculated.

**Figure 11.** The comparison of the FRF of the intact beam and damage structures with one hole and two holes measured by (**a**) piezo-patch MFC1, (**b**) the piezo-patch MFC2.

**Figure 12.** The comparison of power spectrum of the intact beam and damage structures with one hole and two holes measured by (**a**) piezo-patch MFC1, (**b**) the piezo-patch MFC2.

#### **5. Damage Identification Metrics and Discussion**

The numerical and experimental investigations of the intact structure and damage structures presented in the previous sections showed problems with a proper identification of damage in the structure only in terms of FRFs, because the dynamics of these structures is scattered when the frequency increases. Therefore, in order to assess the precision of the damage type and the damage localization, the damage identification metrics should be calculated. Taking this fact into account, five different damage identification metrics M1–M5 taken from [41–45] were calculated for each considered damage structures, and the results of the calculations are presented in this section. In addition, in order to perform better analysis, the metric M2 for the intact structure is calculated as a reference value as well as the damage indicators M1–M5 on the basis of the computational model with three holes. The calculation of damage indices M1, M3, M4, and M5 cannot be done for a healthy baseline because they represent a relationship between the damaged and the intact structure, and consequently a reference M2 marker can link the performances together.

In the first step, their values were calculated based on the frequency responses from the numerical models and then from the FRF generated (Figures 8–10) from the laboratory setup.

The metrices considered in this paper can be divided into two groups: quantitative indicators M1–M3 given by Equations (25)–(27)and qualitative indicators M4–M5 given by Equations (28) and (29). The calculation of these metrics for the first group were performed in terms of a specific frequency value to assess the damage localization in the structure, while in the case of the second group, in terms of the selected frequency range 1–400 Hz, to assess the level of damage. Finally, the results obtained from the computational model considering the damage structure with three holes were collected in Table 3, while results for the experimental response without the structure with the most number of holes are presented in in Table 4. In addition, in order to easier analyze, the damage metric M2 for the intact and damage structures with one and two holes is also presented in the form of a diagram in Figure 13.

$$M\_1 = \max\left( \left. \frac{H\_I(f\_i)}{H\_d(f\_i)} \right|\_{P1}, \left. \frac{H\_I(f\_i)}{H\_d(f\_i)} \right|\_{P2} \right),\tag{25}$$

where:

*HI*(*f*), *Hd*(*f*) denotes frequency response of the intact and damage structure, respectively.

$$M\_2 = \max\left(\frac{H\_d(f\_i)|\_{P2}}{H\_d(f\_i)|\_{P1}}\right),\tag{26}$$

$$M\_3 = \frac{\log\left(H\_D(f\_i) - H\_I(f\_i)\right)}{\log(H\_I(f\_i))} \* 100\%. \tag{27}$$


**Table 4.** Results of the damage metrics M1, M2, and M3 calculated in terms of the generated FRF from the lab stand.



**Table 4.** *Cont.*

**Figure 13.** The comparison of damage identification metric M2 of the damage structures with one hole and two holes calculated based on (**a**) numerical approach, (**b**) experimental approach.

Taking into account the results collected in Tables 3 and 4 and Figure 13, it can be noticed that the experimental tests and the obtained values of the damage metrics M1, M2, and M3 verify in their numerical results. The analysis of the values collected in Tables 3 and 4 show that a gradual decrease of the structure stiffness in a chosen area of the structure leads to the increase of the values of M1 and M3. The inverse effect can be obtained in the case of the analysis of the damage metric M2 (Figure 13), where the increasing number of holes in the structure relative to the intact structure leads to a decrease of its maximum value. It is, however, important to note that all metrics exhibit a monotonicity of calibration against damage, which is important for SHM. Further analysis of these tables also shows that the calculated maximum values of the damage metrics M1 and M3 for the piezo-composite for MFC1 cases are higher than their values calculated for the piezo-patch MFC2. This leads to a conclusion that the damage in the structure is located closer the piezo MFC1 and the piezo MFC2. With adequate sensors, this can lead to the localization of damages better. The actual demand of the spacing of sensors will eventually depend on the demands of detection of the feature of interest in terms of extent and resolution, noise in the measured signal, and the excitation that generates the responses.

Next, the values of the qualitative indicators M4 and M5, given in Equations (28) and (29), were calculated to assess the level of damage in the structure in the selected frequency range of 1–400 Hz, with a frequency increment Δ*f* of 0.00024 Hz (reciprocal of sampling frequency *fs* = 4096 Hz). Similar to previous cases, the first frequency responses functions from the numerical model (Figures 7–9) were taken to calculate these values, and next, the FRF from the laboratory experiment to verify them. The obtained results for the damage metric M4 are shown in Figure 14, while the results for the damage metric M5 in Figure 15.

$$M\_4 = \frac{\Delta f}{f\_{high} - f\_{low}} \sum\_{i=1}^{n} \left| \frac{H\_d(f\_i) - H\_I(f\_i)}{H\_I(f\_i)} \right|,\tag{28}$$

where: Δ*f*—the frequency increment; *fhigh*—the upper frequency; *flow*—the lower frequency.

**Figure 14.** The comparison of damage identification metric M4 of the damage structures with one hole and two holes calculated based on (**a**) numerical approach, (**b**) experimental approach.

**Figure 15.** The comparison of damage identification metric M5 of the damage structures with one hole and two holes calculated based on (**a**) numerical approach, (**b**) experimental approach.

As observed in Figures 14b and 15b, the experimental results of the damage metrics M4 and M5 carried out for only two damage structures with one and two holes allow verifying the values obtained from the computational model. Analysis of M4 indicates a decrease of its values for a gradually increasing number of the drilled holes in a chosen area of the beam even for more damaged structure (see orange line in Figure 14a). A similar observation was noted during the analysis of M5 where its value for the structure with three holes were the lowest. Overall, damage metrics M1–M5 are useful to identify damage and its localization, and can support SHM for beam-like structures. Moreover, these results can be useful to build equivalent damage model and also create a fundamental, low-fidelity system which can lend itself to further studies.

#### **6. Summary and Conclusions**

The use of piezoelectric patches in SHM has expanded the possibilities of use of energy harvesting in recent times. Nonlinearity in the piezo patches with a potential application for SHM has led to investigations in this paper on structures composed of thin piezo-stripes by creating computational models for them. With the current focus on using traditional and modern sensors to aid digital twinning and model updating, such a focus on the behavior of the sensors becomes even more important. Composite structures are making new inroads into a range of sectors, including renewable energy, and so this example is also relevant for future expansion in terms of sustainability of such solutions.

Taking this into account, the stress–strain effect in the laminate structure was analyzed first in this paper to create a fundamental background model. Next, a modal analysis of the chosen structures (intact and damaged) with piezo patch MFCs were carried out using FEM, establishing a homogenized model of MFC elements. Results presented in Figures 4 and 5 indicate that a gradual increase of the number of the holes drilled in the beam in a chosen area slightly affects the values of the resonance and antiresonance frequencies and leads to a decrease of the vibration amplitude due to changes in local stiffness. Further analysis of response to harmonic loading, performed on the FEM models for the intact and damage structures indicates that the distance between the measurement location and excitation leads to collocated (sensor MFC1) and non-collocated (sensor MFC2) systems. Taking the FRF of the collocated system into account (Figures 8a, 9a and 10a) the inversion of resonance and anti-resonance frequencies is observed. For non-collocated system (Figures 8b, 9b and 10b), the increase of distance between the sensor and the location of excitation leads to the estimation of FRFs while omitting the first antiresonance frequency.

Next, FRFs of the FEM models of the intact and damaged structures are used to assess the damage location. Five different damage indexes were calculated in this regard, comprising of three quantitative and two qualitative indicators, estimated as a function of locations of the sensors on the structure. Estimated damage indices in Table 3 show that a gradual decrease of the beam stiffness leads to consistent and monotonic change of the indexes (26% increase of M1, 20–30% increase of M3 and decrease of M2) with respect to the reference value. This consistency of the damage indicators in the presence of realistic conditions is desirable and is indicative of their robustness.

Additionally, taking into account the values in Tables 3 and 4, it can be noticed that the damage of the structure is located closer to the MFC1 sensor and consequently, higher values of damage indicators closer to the harvesting sensors can also be used to identify the approximate location of the damages. With reasonably spaced sensors in the context of damage detection resolution requirements, this can provide information in terms of detection of the damage location. The qualitative metrics also show a decrease in values for increasing number of holes in the structure, as observed from numerical simulations. This consistency of multiple metrics to describe same damage changes also opens up the possibility of using combined metrics to have a more robust detection scheme.

Experimental investigations carried out in the laboratory (one intact and two damage) with an attached PEH allowed to verify the numerical results indicate that estimated FRFs from piezo-sensors MFC1 and MFC2 are consistent with changes due to damage. Subsequent analysis (Figures 11 and 12) also confirms that the beam stiffness decreasing in a chosen region by drilling an increasing number of holes slightly affects the values of resonances and antiresonances, but significantly affects their amplitude within the range of low frequencies. This is especially illustrated in the vicinity of the first natural frequency where drilling subsequent holes leads to the increase of their values (Figure 12a—MFC1), while in the case of the measurement using MFC2 sensor, the amplitude of vibrations increases and then decreases. Heterogenous adhesion between the harvesting elements and the host structure can lead to such a situation. Summarizing the experimental tests, it

can be said that the damage location on the basis only on the analysis of FRF is difficult for energy harvesting and it must be processed further to create relevant markers of damage detection. The proposed damage metrics in this paper illustrate how such markers can be developed and combined, especially in an output-only context. The results of Table 4 indicate similar trends of the proposed metrics as compared to what was observed through numerical simulations. Higher values of damage metrics were observed for sensor closer to the damage (Figures 14b and 15b) along with distinctive and consistent difference over the entire testing range of 1–400 Hz.

This work can act as a reference point for modelling and the expectation of performance of such energy harvesting based SHM sensors for applications in civil/ mechanical systems.

**Author Contributions:** A.K.: concept, supervision, methodology, numerical calculation, experimental investigations, writing of original manuscript; V.P.: writing and editing of the reviewing manuscript; K.L.: determine of computational model for the damage structure with three holes, numerical calculation of damage metrics M1-M5 for the damage structure with three holes. All authors have read and agreed to the published version of the manuscript.

**Funding:** Project financing through the program of the Minister of Science and Higher Education of Poland named "Regional Initiative of Excellence" in 2019–2022 project number 011/RID/2018/19 amount of financing 12,000,000 PLN".

**Institutional Review Board Statement:** Not applicable.

**Informed Consent Statement:** Not applicable.

**Data Availability Statement:** Data for the experiments are available from the authors on request.

**Acknowledgments:** Vikram Pakrashi would like to acknowledge the Science Foundation Ireland MaREI Centre 12/RC/2302\_2 and the UCD Energy Institute. Andrzej Koszewnik would like to acknowledge for participation in the program of the Minister of Science and Higher Education of Poland named "Regional Initiative of Excellence" in 2019–2022 project number 011/RID/2018/19 amount of financing 12,000,000 PLN that is realize on the Faculty of Mechanical Engineering of Bialystok University of Technology.

**Conflicts of Interest:** The authors declare no conflict of interest. The funders had no role in the design of the study; in the collection, analyses, or interpretation of data; in the writing of the manuscript, or in the decision to publish the results.

#### **References**


## *Article* **Experimentally Verified Analytical Models of Piezoelectric Cantilevers in Different Design Configurations**

**Zdenek Machu 1,2,\*, Ondrej Rubes 1, Oldrich Sevecek <sup>1</sup> and Zdenek Hadas <sup>1</sup>**


**Abstract:** This paper deals with analytical modelling of piezoelectric energy harvesting systems for generating useful electricity from ambient vibrations and comparing the usefulness of materials commonly used in designing such harvesters for energy harvesting applications. The kinetic energy harvesters have the potential to be used as an autonomous source of energy for wireless applications. Here in this paper, the considered energy harvesting device is designed as a piezoelectric cantilever beam with different piezoelectric materials in both bimorph and unimorph configurations. For both these configurations a single degree-of-freedom model of a kinematically excited cantilever with a full and partial electrode length respecting the dimensions of added tip mass is derived. The analytical model is based on Euler-Bernoulli beam theory and its output is successfully verified with available experimental results of piezoelectric energy harvesters in three different configurations. The electrical output of the derived model for the three different materials (PZT-5A, PZZN-PLZT and PVDF) and design configurations is in accordance with lab measurements which are presented in the paper. Therefore, this model can be used for predicting the amount of harvested power in a particular vibratory environment. Finally, the derived analytical model was used to compare the energy harvesting effectiveness of the three considered materials for both simple harmonic excitation and random vibrations of the corresponding harvesters. The comparison revealed that both PZT-5A and PZZN-PLZT are an excellent choice for energy harvesting purposes thanks to high electrical power output, whereas PVDF should be used only for sensing applications due to low harvested electrical power output.

**Keywords:** energy harvesting; vibrations; piezoelectric; analytical model; beam model; equivalent model; power prediction

#### **1. Introduction**

Energy harvesting is more than 20 years a hot topic in the field of wireless sensing [1] since it allows for converting various energy types from ambient sources into an electrical one. Although the amount of such harvested energy is usually small (tens of μW up to several mW), it can be used as a source of electrical power for modern, low power-consuming sensors that are typically used in wearable electronics and industrial applications [2] where powering using cables is not feasible (either due to a hazardous environment or complex setup). Piezoelectric kinetic energy harvesters in the form of a vibrating multilayer structure with piezoelectric layers [3] are commonly used in vibration energy harvesting applications, where the structure is excited by an ambient source of vibrations. The main task of kinetic energy harvesters is then to transform the mechanical energy of ambient vibrations, mainly those of machine frames or human body movement, into useful electrical energy by means of the direct piezoelectric phenomenon.

The main goal in the field of energy harvesting is to design a kinetic energy harvester which is capable to generate a sufficient amount of electrical energy in a particular vibratory

**Citation:** Machu, Z.; Rubes, O.; Sevecek, O.; Hadas, Z. Experimentally Verified Analytical Models of Piezoelectric Cantilevers in Different Design Configurations. *Sensors* **2021**, *21*, 6759. https:// doi.org/10.3390/s21206759

Academic Editor: Salvatore Salamone

Received: 17 September 2021 Accepted: 9 October 2021 Published: 12 October 2021

**Publisher's Note:** MDPI stays neutral with regard to jurisdictional claims in published maps and institutional affiliations.

**Copyright:** © 2021 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (https:// creativecommons.org/licenses/by/ 4.0/).

environment [4] in order to power some other electronic equipment. However, each application has its requirements or limits for dimensions and weight of the harvester; the principle of energy harvesting can be used practically everywhere, for example in the field of medicine [5], wearables [6], portables [7], aircrafts [8], structural health monitoring of railways [9] or bridges [10].

It has been proved many times that for harvesting energy from ambient vibrations the kinematically excited cantilever beam is one of the most effective designs of a piezoelectric energy harvester. The fundamental and also the most important issue of this solution is the choice of a suitable piezoelectric material for effective electromechanical conversion. The review of commonly used piezoelectric materials and structures for energy harvesting purposes is summarized in publication [11], where it is shown that not only the material itself but also the intended operational mode significantly affects the amount of harvested power due to a great variation in piezoelectric coefficients. The highest piezoelectric coefficients (generally, the higher the coefficients, the higher the amount of harvested power) are provided by piezoceramic materials [12], especially those based on lead (PZT). As a non-toxic alternative, new lead-free piezoceramic materials have been developed which are based on multifunctional Perovskite [13] or structured layers made of Barium and Titanate [14]. Besides these piezoceramic materials which are inherently very brittle and stiff, there are also more flexible materials such as macro-fiber composites which are very promising in the area of strain energy harvesting [15] and piezopolymers which are summarized in review paper [16]. An example of a cantilever harvesting device based on a piezoelectric polymer (PVDF) is presented in paper [17] and the effectivity of PVDF in energy harvesting applications is nowadays widely discussed [6].

In conclusion, the two most important factors that determine the effectiveness of a vibrational energy harvesting device are the used piezoelectric material and the harvester's geometry. Many recent works were concerned about the optimal harvester's geometry for selected piezoelectric material, e.g., [18], but the effectivity of various piezoelectric materials has not been widely discussed yet. Both the selection of efficient piezoelectric material and suitable geometry of the harvester can be solved with an appropriate model of the piezoelectric resonator. For this reason, the presented paper is organized as follows. First, derivation of an analytical beam model of a kinematically excited piezoelectric cantilever in both bimorph/unimorph configurations which also respects the dimensions of used tip mass. This beam model is subsequently reduced to a single degree-of-freedom (DOF) system using the first mode shape function. Although, the derivation of a coupled electromechanical model was published several times, e.g., [19–23], here, we also show the effect of chosen mode shape function which is used in reducing the beam model into single DOF model. Then, the model is verified with 3 different experimental results. Finally, the main aim of this paper is to provide a methodology based on a verified model that can be used to compare the effectivity of materials commonly used in energy harvesting applications.

#### **2. Model of Piezoelectric Vibration Energy Harvester**

In order to harvest as much energy from vibrations as possible, it is paramount to properly design dimensions of the harvester and optimize its electrical impedance. This goal can easily be achieved with an analytical model which is able to predict electromechanical response of piezoelectric energy harvesters. Therefore, here in this paper a single DOF model of a kinematically excited cantilever in both bimorph/unimorph configurations is derived. This analytical model is based on Euler-Bernoulli (thin) beam theory and its output is compared with results from experiments conducted with three different piezoelectric harvesters described further in the following section. Then, the derived model is used in a comparative study to compare the piezoelectric materials used in the experiments in terms of energy harvesting efficiency.

#### *2.1. Bimorph Cantilever Beam with Piezoelectric Layers in Series*

A piezoelectric cantilever beam harvester in a bimorph configuration is shown in Figure 1. This beam model with dimensions *L* × *B* × *H*, where *H* = 2 *h*<sup>p</sup> + *h*s, is used for obtaining a single DOF analytical model. The clamped end is kinematically excited with a time-harmonic base acceleration *a*(*t*) from an external source of vibrations. Piezoelectric layers are in operational mode 31 (in-thickness polarization of piezoelectric layers and axial bending deformation of the harvester) whose polarization is denoted by arrow symbols in the figure. These layers have electrodes present over a region of dimensions *L*<sup>E</sup> × *B* which is mentioned further in the text as section *V*E; the remaining portion of piezoelectric layers is not polarized (mentioned as section *V*R), and thus is not affected by the piezoelectric effect. A tip mass *M*<sup>t</sup> of negligible rotary inertia is attached to the free end of the beam spanning over the length *L*Mt—this section is denoted as *V*R, Mt. The bimorph model is reduced to a single DOF model which describes the movement of the bimorph's free end *q* relative to the moving clamped end.

**Figure 1.** Geometric model of a piezoelectric bimorph in operational mode 31.

The assumptions of Euler-Bernoulli beam theory combined with those of the classical laminate theory concerning continuous strain throughout the layers of a multilayer structure imply that the strain *ε<sup>x</sup>* within bimorph's layers can be expressed as

$$
\varepsilon\_x = -\frac{\partial^2 w}{\partial x^2} z\_\prime \tag{1}
$$

where *w* is transverse displacement of the beam's centerline relative to the movement of the excited clamped end. Since the beam is split into a section *V*<sup>E</sup> which is affected by the piezoelectric effect, section *V*<sup>R</sup> which is not affected by the piezoelectric effect and section *V*R, Mt with distributed tip mass attached, the displacement needs to be a piecewise function defined as

$$w(\mathbf{x},t) = \begin{cases} \begin{array}{c} w\_1(\mathbf{x},t) \text{ for } \mathbf{x} \in [0,L\_E] \\ w\_2(\mathbf{x},t) \text{ for } \mathbf{x} \in (L\_{E'},L - L\_{Mt}] \\ w\_3(\mathbf{x},t) \text{ for } \mathbf{x} \in (L\_- - L\_{Mt},L] \end{array} . \end{cases} . \tag{2}$$

For these functions it must hold that they are continuous and smooth at intersection points, i.e., they must satisfy following conditions:

$$\begin{aligned} w\_1(\mathbf{x} = L\_{E\prime} \ t) &= w\_2(\mathbf{x} = L\_{E\prime} \ t) \\ w\_2(\mathbf{x} = L - L\_{\mathbf{M}\prime} \ t) &= w\_3(\mathbf{x} = L - L\_{\mathbf{M}\prime} \ t) \\ \frac{\partial w\_1}{\partial \mathbf{x}}(\mathbf{x} = L\_{E\prime} \ t) &= \frac{\partial w\_2}{\partial \mathbf{x}}(\mathbf{x} = L\_{E\prime} \ t) \\ \frac{\partial w\_2}{\partial \mathbf{x}}(\mathbf{x} = L - L\_{\mathbf{M}\prime} \ t) &= \frac{\partial w\_3}{\partial \mathbf{x}}(\mathbf{x} = L - L\_{\mathbf{M}\prime} \ t) \end{aligned} \tag{3}$$

The piezoelectric constitutive relations, which are required further in the derivation of the analytical model, take the following form for the uniaxial stress state [23]:

$$
\varepsilon\_{\mathbf{x}} = \mathbf{s}\_{11}^{\mathrm{E}} \sigma\_{\mathbf{x}} + d\_{\Im 1} E\_{\mathbf{z} \prime} \tag{4}
$$

$$D\_z = d\_{31} \sigma\_x + \epsilon\_{33}^T E\_{z\prime} \tag{5}$$

where *ε<sup>x</sup>* represents normal strain, *s<sup>E</sup>* <sup>11</sup> is mechanical compliance measured at constant electric field, *σ<sup>x</sup>* is normal stress, *d*<sup>31</sup> is the mode-31 component of the piezoelectric charge coefficient matrix, *Ez* is the *z*-component of electric field intensity, *Dz* is the *z*-component of electric flux density and *<sup>T</sup>* <sup>33</sup> is permittivity in direction of the polarization axis measured at constant mechanical stress. From Equation (4) the expression for the stress *σ<sup>x</sup>* is extracted as

$$
\sigma\_x = \left(s\_{11}^E\right)^{-1} \varepsilon\_x - \left(s\_{11}^E\right)^{-1} d\_{31} E\_z \,. \tag{6}
$$

Note that the reciprocal value of *s<sup>E</sup>* <sup>11</sup> equals to the elastic modulus *Yp* of used piezoelectric material. Equation (6) can then be rewritten as

$$
\sigma\_x = \mathcal{Y}\_p \varepsilon\_x - \underbrace{\mathcal{Y}\_p d\_{31}}\_{\varepsilon\_{31}} E\_z = \mathcal{Y}\_p \varepsilon\_x - \varepsilon\_{31} E\_z \,\,\,\,\tag{7}
$$

where *e*<sup>31</sup> is the piezoelectric modulus. Substituting (7) into (5) yields

$$D\_z = \underbrace{d\_{31} Y\_p \varepsilon\_x + \underbrace{\left(\varepsilon\_{33}^T - d\_{31} \varepsilon\_{31}\right)}\_{\varepsilon\_{33}^S} E\_z = \varepsilon\_{31} \varepsilon\_x + \epsilon\_{33}^S E\_{z\prime} \tag{8}$$

where *<sup>S</sup>* <sup>33</sup> is the permittivity of used piezoelectric material measured at constant strain.

Piezoelectric materials are dielectrics and, as a consequence of the Gauss' law, it holds that *∂Dz*/*∂z* =0[24]. This implies that *Dz* = const. Since the strain term in (8) changes linearly with the *z*-coordinate, we can use its mean value at the center of the *n*-th piezoelectric layer; the *z*-coordinate of the center of *n*-th layer is denoted as *z*Tp*n*. Next, the fundamentals of electricity, see e.g., [25], state that integration of *Ez* over the thickness of a particular layer yields the voltage drop for the given layer. In order to make *Dz* independent of the *z*-coordinate, *Ez* must be a linear function of this coordinate with its mean value at the center of *n*-th layer given by

$$E\_z \left( z\_{\rm T^u} \right) = -\frac{\mathcal{U}\_n}{h\_p} = -\frac{\mathcal{U}}{2h\_p} \tag{9}$$

where *U* is the magnitude of generated voltage drop and *hp* is the thickness of piezoelectric layers. By using Equations (1) and (9) in Equation (8) and by assuming that *Dz* is a layer-wise function, one receives:

$$D\_z = -\epsilon\_{31} z\_{\rm Tp} \frac{\partial^2 w}{\partial x^2} - \frac{\epsilon\_{33}^S}{2h\_{\rm P}} \mathcal{U}. \tag{10}$$

By combining Equations (8) and (10), the following expression for *Ez* is obtained:

$$E\_z = \frac{c\_{31}}{\epsilon\_{33}^S} \frac{\partial^2 w}{\partial \mathbf{x}^2} \left( z - z\_{\rm TP^n} \right) - \frac{\mathcal{U}}{2h\_p}. \tag{11}$$

According to [26], the electric current *I* generated by two in-series connected piezoelectric layers can be expressed as

$$I = \frac{\mathcal{U}}{\mathcal{R}\_l} = \frac{\mathcal{d}}{\mathcal{d}t} \left( \iint\_{A\_E} D\_z \mathbf{dx} dy \right), \tag{12}$$

where the subscript *A*<sup>E</sup> denotes the area of electrodes and *R*<sup>l</sup> represents the connected resistance. Inserting (10) into Equation (12) yields

$$\frac{dI}{d\mathbf{k}\_{\parallel}} = \frac{\mathbf{d}}{\mathbf{d}t} \left[ \iint\_{A\_{E}} \left( -\varepsilon\_{31} z\_{\Gamma \text{Pl}} \frac{\partial^{2} w}{\partial \mathbf{x}^{2}} - \frac{\varepsilon\_{33}^{S}}{2h\_{\text{P}}} \mathbf{d}I \right) \mathbf{dx} \mathbf{dy} \right]. \tag{13}$$

Integrating terms in (13) leads to a following PDE which governs the electrical behavior of the considered bimorph

$$\mathbf{C}\_{\rm eq} \frac{\mathbf{d} \mathbf{U}}{\mathbf{d}t} + \frac{1}{R\_l} \mathbf{U} = \kappa \frac{\partial^2 w}{\partial \mathbf{x} \partial t} \Big|\_{\mathbf{x} = \mathbf{L}\_{\rm E}} \tag{14}$$

where *C*eq denotes the bimorph's equivalent capacitance defined as

$$\mathcal{L}\_{\text{eq}} = \epsilon\_{33}^{\text{S}} \frac{BL\_{\text{E}}}{2h\_{\text{P}}},\tag{15}$$

and *κ* represents generic electromechanical coupling defined as

$$
\kappa = -B\varepsilon\_{31} z\_{\rm T\_{\rm P}}.\tag{16}
$$

As a next step, Equations (11) and (1) can be inserted into Equation (7) to obtain an expression for stress *σ<sup>x</sup>* within the polarized piezoelectric layers:

$$
\sigma\_x = -\mathcal{Y}\_\mathsf{P} \frac{\partial^2 w}{\partial \mathbf{x}^2} z + \frac{e\_{31}^2}{\epsilon\_{33}^S} \frac{\partial^2 w}{\partial \mathbf{x}^2} (z\_{\mathsf{T}\mathsf{P}^n} - z) + \frac{e\_{31}}{2h\_\mathsf{P}} \mathsf{U}. \tag{17}
$$

Within the other layers (non-polarized piezoelectric layers and the substrate) the stress obeys the Hooke's law:

$$
\sigma\_{\mathbf{x}} = -Y\_{\mathbf{n}} \frac{\partial^2 w}{\partial \mathbf{x}^2} \mathbf{z} \,, \tag{18}
$$

where *Y*n denotes elastic modulus of the used piezoelectric material *Y*p or the substrate *Y*s. Total energy stored in the considered bimorph upon vibrations consists of kinetic energy *E*k, strain energy *E*<sup>p</sup> and the work done by inertial forces due to kinematic excitation *W*ext. Kinetic energy of the considered beam can be written as

$$\begin{split} E\_{\mathbf{k}} &= \frac{1}{2} \iiint\_{\mathrm{V}} \left[ \rho\_{\mathrm{n}} \left( \frac{\partial w\_{\mathrm{l}}}{\partial t} \right)^{2} \right] \mathrm{d}V + \frac{1}{2} \iiint\_{\mathrm{V}} \left[ \rho\_{\mathrm{n}} \left( \frac{\partial w\_{\mathrm{2}}}{\partial t} \right)^{2} \right] \mathrm{d}V + \frac{1}{2} \iiint\_{\mathrm{V}\_{\mathrm{sh}}} \left[ \left( \rho\_{\mathrm{n}} + \frac{\mathrm{M}\_{\mathrm{t}}}{\mathrm{L} \mathrm{M} \mathrm{H}} \right) \left( \frac{\partial w\_{\mathrm{3}}}{\partial t} \right)^{2} \right] \mathrm{d}V \\ &= \frac{1}{2} \iiint\_{\mathrm{V}} \left[ m^{\*} \left( \frac{\partial w\_{\mathrm{1}}}{\partial t} \right)^{2} \right] \mathrm{d}x + \frac{1}{2} \iiint\_{\mathrm{L}\_{\mathrm{E}}} \left[ m^{\*} \left( \frac{\partial w\_{\mathrm{2}}}{\partial t} \right)^{2} \right] \mathrm{d}x + \frac{1}{2} \iiint\_{\mathrm{L}\_{\mathrm{t}}} \left[ \left( m^{\*} + \frac{\mathrm{M}\_{\mathrm{t}}}{\mathrm{L}\_{\mathrm{hh}}} \right) \left( \frac{\partial w\_{\mathrm{3}}}{\partial t} \right)^{2} \right] \mathrm{d}x, \end{split} \tag{19}$$

where *ρ<sup>n</sup>* is density of the *n*-th layer, *M*<sup>t</sup> is the attached tip mass and *m*\* is the bimorph's mass per unit of its length defined as

$$m^\* = B\left(\rho\_{\mathfrak{s}} h\_{\mathfrak{s}} + 2\rho\_{\mathfrak{p}} h\_{\mathfrak{p}}\right). \tag{20}$$

Strain energy stored in the bimorph can be expressed as

$$E\_{\rm P} = \frac{1}{2} \iiint\_{\rm E} [\varepsilon\_{\rm F} \sigma\_{\rm x}] \, \mathrm{d}V + \frac{1}{2} \iiint\_{\rm V\_{\rm R}} \left[ Y\_{\rm R} \varepsilon\_{\rm x}^{2} \right] \, \mathrm{d}V + \frac{1}{2} \iiint\_{\rm V\_{\rm R}} \left[ Y\_{\rm R} \varepsilon\_{\rm x}^{2} \right] \, \mathrm{d}V$$

$$= \frac{1}{2} \int\_{0}^{L\_{\rm E}} \left[ I\_{\rm Pdxo}^{\*} \left( \frac{\partial^{2} w\_{\rm i}}{\partial x^{2}} \right)^{2} - \kappa w\_{1} \frac{q\mathcal{L}}{\mathrm{d}x} (\mathbf{x} - L\_{\rm E}) \mathrm{d}l \right] \mathrm{d}\mathbf{x} + \frac{1}{2} \int\_{L\_{\rm E}}^{L - L\_{\rm M}} \left[ I^{\*} \left( \frac{\partial^{2} w\_{\rm i}}{\partial x^{2}} \right)^{2} \right] \mathrm{d}\mathbf{x} + \frac{1}{2} \int\_{L - \frac{1}{2} \mathrm{d}l \mathrm{d}t}^{\dagger} \left[ I^{\*} \left( \frac{\partial^{2} w\_{\rm i}}{\partial x^{2}} \right)^{2} \right] \mathrm{d}\mathbf{x} \right. \tag{21}$$

where dδ/d*x* is the first derivative of Dirac's delta function, *J*∗ piezo denotes bending stiffness of the beam section where the polarization of a piezoelectric material is considered (over the length *L*E) defined as

$$J\_{\rm piezo}^{\star} = \frac{1}{12} \mathcal{Y}\_{\rm s} B h\_{\rm s}^{3} + 2 \mathcal{Y}\_{\rm P} \left[ \frac{1}{12} B h\_{\rm P}^{3} + \left( \frac{h\_{\rm P}}{2} + \frac{h\_{\rm s}}{2} \right)^{2} B h\_{\rm P} \right] + \upsilon\_{\rm } \tag{22}$$

where

$$\upsilon = 2 \times \frac{1}{12} \frac{e\_{31}^2}{e\_{33}^8} Bh\_{\mathbb{P}}^3 \, . \tag{23}$$

and *J* \* is bending stiffness of the non-polarized section of the beam (the rest of the beam outside the length *L*E) defined as

$$J^\* = \frac{1}{12} Y\_\text{s} B h\_\text{s}^3 + 2Y\_\text{P} \left[ \frac{1}{12} B h\_\text{P}^3 + \left( \frac{h\_\text{P}}{2} + \frac{h\_\text{s}}{2} \right)^2 B h\_\text{P} \right]. \tag{24}$$

The work done by inertia forces due to kinematic excitation is defined as

$$\begin{array}{rcl} \mathsf{W\_{ext}} & = -\iiint\_{\mathrm{E}} [\rho\_{n}a\_{0}w\_{1}] \mathrm{d}V - \underset{\mathrm{V\_{R}}}{\mathrm{III}} [\rho\_{n}a\_{0}w\_{2}] \mathrm{d}V - \underset{\mathrm{V\_{R,\mathrm{Mt}}}}{\mathrm{III}} \left[ \left( \rho\_{n} + \frac{M\_{\mathrm{t}}}{L\_{\mathrm{Mt}}\mathrm{H}t} \right) a\_{0}w\_{3} \right] \mathrm{d}V \\ & = -\underset{\mathrm{D}}{\mathrm{I}} [m^{\*}a\_{0}w\_{1}] \mathrm{d}\mathbf{x} - \underset{\mathrm{L\_{E}}}{\mathrm{I}} [m^{\*}a\_{0}w\_{2}] \mathrm{d}\mathbf{x} - \underset{\mathrm{L\_{-}}}{\int\_{\mathrm{L}}} \left[ \left( m^{\*} + \frac{M\_{\mathrm{t}}}{L\_{\mathrm{Mt}}} \right) a\_{0}w\_{3} \right] \mathrm{d}\mathbf{x}. \end{array} \tag{25}$$

Subsequently, Hamilton's variational principle [27] is used to obtain equations of motion in the form of PDEs with a nonzero right-hand side. The equations of motion for the polarized portion of the beam and for the non-polarized portions of the beam take the following form:

$$J\_{\rm plezo}^{\*} \frac{\partial^{4} w\_{1}}{\partial \mathbf{x}^{4}} + m^{\*} \frac{\partial^{2} w\_{1}}{\partial t^{2}} - \kappa \frac{\mathbf{d} \delta}{\mathbf{dx}} (\mathbf{x} - L\_{\rm E}) \mathcal{U} = -m^{\*} a\_{0 \prime} \ge c \left[ 0, \ L\_{\rm E} \right] \tag{26}$$

$$J^\* \frac{\partial^4 w\_2}{\partial x^4} + m^\* \frac{\partial^2 w\_2}{\partial t^2} = -m^\* a\_{0\prime} \ge \varepsilon (L\_{\rm E\prime} \ L - L\_{\rm M\rm t}) \tag{27}$$

$$J^\* \frac{\partial^4 w\_3}{\partial x^4} + \left(m^\* + \frac{M\_\text{t}}{L\_\text{Mt}}\right) \frac{\partial^2 w\_3}{\partial t^2} = -\left(m^\* + \frac{M\_\text{t}}{L\_\text{Mt}}\right) a\_{0\prime} \ge \in (L - L\_\text{Mt}, L] \tag{28}$$

Equations above, however, do not account for damping; therefore, they have to be extended with a damping term. Here, we shall consider the stiffness damping term from Rayleigh's Damping theorem:

$$J\_{\rm piezo}^{\*} \frac{\partial^4 w\_1}{\partial \mathbf{x}^4} + \frac{2b\_r}{\Omega\_1} l\_{\rm piezo}^{\*} \frac{\partial^5 w\_1}{\partial \mathbf{x}^4 \partial t} + m^\* \frac{\partial^2 w\_1}{\partial t^2} - \mathbf{x} \frac{\mathbf{d}\delta}{\mathbf{dx}} (\mathbf{x} - L\_\mathbf{E}) \mathcal{U} = -m^\* a\_{0\prime} \ge \in [0, L\_\mathcal{E}] \tag{29}$$

$$J^\* \frac{\partial^4 w\_2}{\partial \mathbf{x}^4} + \frac{2b\_r}{\Omega\_1} J^\* \frac{\partial^5 w\_2}{\partial \mathbf{x}^4 \partial t} + m^\* \frac{\partial^2 w\_2}{\partial t^2} = -m^\* a\_{0\prime} \ge \in (L\_\text{E}, L \ge -L\_{\text{Mft}}] \tag{30}$$

$$J^\* \frac{\partial^4 w\_3}{\partial x^4} + \frac{2b\_\mathbf{r}}{\Omega\_1} l^\* \frac{\partial^5 w\_3}{\partial x^4 \partial t} + \left(m^\* + \frac{M\_\mathbf{l}}{L\_{\text{Mft}}}\right) \frac{\partial^2 w\_3}{\partial t^2} = -\left(m^\* + \frac{M\_\mathbf{l}}{L\_{\text{Mft}}}\right) a\_{0\prime} \ge \in (L - L\_{\text{Mft}}, L] \tag{31}$$

where *b*<sup>r</sup> is the considered damping ratio and Ω<sup>1</sup> is the value of the beam's first eigenfrequency. Equations (29)–(31) together with (14) form a complete equation system which describes the electromechanical response of the considered bimorph.

However, this system of PDEs is actually not very effective to be used in modelling of energy harvesting devices because of its complexity, thus its transformation into a much simpler single DOF model is necessary. In the scope of vibrational energy harvesting applications, the beam is kinematically excited with frequencies very close or equal to the harvester's first resonant frequency *f* 1,r. This fact means that beam vibrations are composed mostly of the first vibrational mode and, as a consequence, the beam's displacement relative to the base movement in all sections (*V*E, *V*<sup>R</sup> and *V*R, Mt) can be written as

$$w(x,t) \approx \phi\_1(x)\eta\_1(t),\tag{32}$$

where *φ*<sup>1</sup> is the mode shape function of the first mode and *η*<sup>1</sup> is its modal coordinate. The shape function *φ*<sup>1</sup> can be approximated with an arbitrary function that resembles the shape of the first bending mode. Although Erturk in [26] recommends using an approximative function which accounts for a tip mass at the beam's free end, such a function is not appropriate for tip masses spanning over a finite length of the beam. Therefore, the following expression was chosen to simplify his approximative function into:

$$\phi\_1(\mathbf{x}) = \mathcal{C}\_1 \cdot \left[ \cos \frac{\lambda\_1}{L} \mathbf{x} \; - \; \cos \mathbf{h} \frac{\lambda\_1}{L} \mathbf{x} + \zeta\_1 \left( \sin \frac{\lambda\_1}{L} \mathbf{x} \; - \; \sin \mathbf{h} \frac{\lambda\_1}{L} \mathbf{x} \right) \right],\tag{33}$$

where

$$\zeta\_1 = \frac{\sin \lambda\_1 - \sin \mathbf{h} \lambda\_1}{\cos \lambda\_1 + \cos \mathbf{h} \lambda\_1}. \tag{34}$$

The eigenvalue *λ*<sup>1</sup> is obtained as the first positive root of the following transcendental equation [26]

$$1 + \cos \lambda\_1 \cos h \lambda\_1 = 0.\tag{35}$$

Equation (33) accurately describes the first mode shape of a beam without a tip mass. Further in the paper it will be shown that simpler functions which deviate from the actual shape overestimate the beam's stiffness and influence the calculated results.

The constant *C*<sup>1</sup> in (33) should be evaluated so that *φ*<sup>1</sup> is mass-normalized to prevent numerical errors in further calculations of the model's parameters, i.e., *φ*<sup>1</sup> satisfies the following condition

$$\int\_{0}^{L-L\_{\text{Mt}}} \phi\_{1}(\mathbf{x}) m^{\*} \phi\_{1}(\mathbf{x}) \mathbf{dx} + \int\_{L-L\_{\text{Mt}}}^{L} \phi\_{1}(\mathbf{x}) \left(m^{\*} + \frac{M\_{\text{t}}}{L\_{\text{Mt}}}\right) \phi\_{1}(\mathbf{x}) \mathbf{dx} = 1 \tag{36}$$

Then, approximation (32) can be inserted into the equation system (14), (29)–(31) which can now be solved effectively using the Galerkin method [27], resulting into a much simpler equation system:

$$M\frac{\mathbf{d}^2\eta}{\mathbf{d}t^2} + B\frac{\mathbf{d}\eta}{\mathbf{d}t} + K\eta + \theta \mathbf{U} = F\_\prime \tag{37}$$

$$\mathbf{C}\_{\rm eq} \frac{\mathbf{d} \mathbf{U}}{\mathbf{d}t} + \frac{1}{R\_I} \mathbf{U} = \theta \frac{\mathbf{d} \eta}{\mathbf{d}t} \,\tag{38}$$

where

$$\begin{aligned} \mathcal{M} &= \int\_0^{L-L\_{\text{Mth}}} \Phi\_1(\mathbf{x}) m^\* \Phi\_1(\mathbf{x}) \mathbf{dx} + \int\_{L-L\_{\text{Mth}}}^L \Phi\_1(\mathbf{x}) \left( m^\* + \frac{M\_1}{L\_{\text{Mth}}} \right) \Phi\_1(\mathbf{x}) \mathbf{dx} = 1 \\ K &= J\_{\text{piezoo}}^\* \int\_0^{L\_{\text{p}}} \left[ \frac{\mathbf{d}^2 \phi\_1(\mathbf{x})}{\mathbf{dx}^2} \right]^2 \mathbf{dx} + J^\* \int\_{L\_{\text{p}}}^L \left[ \frac{\mathbf{d}^2 \phi\_1(\mathbf{x})}{\mathbf{dx}^2} \right]^2 \mathbf{dx} = \Omega\_1^2 \end{aligned} \tag{39}$$
 
$$B = 2h\_r \Omega\_1$$
 
$$\theta = \mathbf{x} \frac{\mathbf{d} \phi\_1}{\mathbf{dx}} \Big|\_{\mathbf{x} = L\_{\text{E}}} $$
 
$$F = - \begin{bmatrix} m^\* a\_0 \int\_0^{L-L\_{\text{Mth}}} \Phi\_1(\mathbf{x}) \mathbf{dx} + \left( m^\* + \frac{M\_1}{L\_{\text{Mth}}} \right) a\_0 \int\_{L-L\_{\text{Mth}}}^L \Phi\_1(\mathbf{x}) \mathbf{dx} \end{bmatrix}$$

#### 2.1.1. Effect of Chosen Mode Shape Function on Model Output

This section addresses the effect of the chosen approximative mode shape function *φ*<sup>1</sup> on the model's behavior. To this purpose, a reference configuration of a piezoelectric harvester with a significant tip mass is needed. This requirement is satisfied by a PZT-5A bimorph from a well-known work of Erturk and Inman [26].

To analyze the influence of the chosen approximative function *φ*<sup>1</sup> on the model's output, the actual mode shape of the reference harvester is needed. To obtain the actual mode shape, a 3D numerical model of the reference harvester was created in commercial FE software ANSYS APDL made of approx. 1,000 SOLID186 higher-order elements. The actual mode shape denoted as *φ*1,true was obtained from a modal analysis of the model using the path post-processing tool. Both the actual mode shape *φ*1,true and the approximation *φ*1,approx defined by (33) normed to unity are plotted in Figure 2a. While mode shapes in the graph look almost identical, a much clearer distinction can be seen in Figure 2b by plotting their first derivatives with respect to *x* (the slope of the mode shape). Here, the actual mode shape *φ*1,true shows a much higher degree of compliance (higher value of d*φ*1/d*x*) at the beam's free end, which is crucial for high power output. This increase in compliance at the beam's free end is caused by the presence of the tip mass. Therefore, the presence of heavy tip masses at the beam's free end causes an increase in beam's compliance near its free end which cannot be accounted for using simpler approximative functions, such as polynomials. Using simpler approximative functions will lead to stiffer behavior of the beam model and result in higher resonant frequencies and underestimation of generated electrical power. Nevertheless, the errors in the model's output by using (33) are not significant as demonstrated further in the paper.

**Figure 2.** (**a**) Comparison between approximation and true mode shape; (**b**) comparison of slopes between approximation and true mode shape.

#### 2.1.2. Single DOF Model of Bimorph Configuration

The system of Equations (37) and (38) is still not suitable for prediction of harvested power due to using the modal coordinate. For this reason, these equations are transformed from the modal coordinate into direct calculation of the relative movement of the bimorph's free end *q* defined by (32) as

$$q(t) = \phi\_1(\mathbf{x} = L)\eta\_1(t),\tag{40}$$

Inserting (40) into (37) and (38) transforms the equation system into the sought single DOF model:

$$M\_{\rm eff} \frac{\mathrm{d}^2 q}{\mathrm{d}t^2} + B\_{\rm eff} \frac{\mathrm{d}q}{\mathrm{d}t} + K\_{\rm eff} q + \theta\_{\rm eff} lI = F\_{\rm eff} \tag{41}$$

$$\text{C}\_{\text{eq}} \frac{\text{d} \text{d} \text{I}}{\text{d} \text{t}} + \frac{1}{R\_{\text{I}}} \text{d} \text{I} = \left. \theta\_{\text{eff}} \frac{\text{d} \text{q}}{\text{d} \text{t}} \right| \tag{42}$$

where

$$\begin{aligned} M\_{\mathrm{eff}} &= \frac{M}{\Phi\_1^2(\mathbf{x} = L)} = \frac{1}{\Phi\_1^2(\mathbf{x} = L)}\\ K\_{\mathrm{eff}} &= \frac{K}{\phi\_1^2(\mathbf{x} = L)} = \frac{\Omega\_1^2}{\phi\_1^2(\mathbf{x} = L)}\\ B\_{\mathrm{eff}} &= \frac{B}{\phi\_1^2(\mathbf{x} = L)} = \frac{2b\_1 \Omega\_1}{\phi\_1^2(\mathbf{x} = L)}\\ \theta\_{\mathrm{eff}} &= \frac{\theta}{\phi\_1(\mathbf{x} = L)} = \frac{\kappa}{\phi\_1(\mathbf{x} = L)} \left. \frac{d\Phi\_1}{\mathrm{d\mathbf{x}}} \right|\_{\mathbf{x} = L\_{\mathrm{E}}}\\ F\_{\mathrm{eff}} &= \frac{F}{\phi\_1(\mathbf{x} = L)} \end{aligned} \tag{43}$$

where capacitance *C*eq is defined through Equation (15) and the connected resistive load *R*<sup>l</sup> represents the useful electrical load.

#### *2.2. Modification of Single DOF Model for Unimorph Configuration*

The model of a unimorph configuration of a piezoelectric harvester which considers only one piezoelectric layer is commonly used with piezoelectric polymers. The unimorph geometric model shares the same parameters to that of a bimorph shown in Figure 1. The single DOF model for the unimorph uses exactly the same equations that were derived for the bimorph, i.e., (41) and (42). Contrary to the bimorph, however, the major difference lies in fact that the bimorph's neutral axis is coincident with its geometrical midplane, whereas this is not true in case of a unimorph. Therefore, the coefficients in (41) and (42) have to be re-defined to respect this fact.

First, the neutral axis of the unimorph *z*<sup>N</sup> is calculated as [28]

$$z\_{\rm N} = \frac{Y\_{\rm s} z\_{\rm Ts} h\_{\rm s} + Y\_{\rm P} z\_{\rm TP} h\_{\rm P}}{Y\_{\rm s} h\_{\rm s} + Y\_{\rm P} h\_{\rm P}} = \frac{1}{2} h\_{\rm s} h\_{\rm P} \frac{Y\_{\rm P} - Y\_{\rm s}}{Y\_{\rm s} h\_{\rm s} + Y\_{\rm P} h\_{\rm P}}.\tag{44}$$

Then, the coefficients in (41) and (42) are re-calculated with respect to the unimorph's neutral axis *z*<sup>N</sup> using the same shape function *φ*<sup>1</sup> as in (33). First, the mass coefficient *M*eff is defined as

$$M\_{\rm eff} = \frac{\int\_0^{L - L\_{\rm Mt}} \phi\_1(\mathbf{x}) m^\* \phi\_1(\mathbf{x}) d\mathbf{x} + \int\_{L - L\_{\rm Mt}}^L \phi\_1(\mathbf{x}) \left(m^\* + \frac{M\_t}{L\_{\rm Mt}}\right) \phi\_1(\mathbf{x}) d\mathbf{x}}{\phi\_1^2(\mathbf{x} = L)},\tag{45}$$

where *m*\* changes to

$$m^\* = B\left(\rho\_{\mathfrak{s}} h\_{\mathfrak{s}} + \rho\_{\mathbb{P}} h\_{\mathbb{P}}\right). \tag{46}$$

Then, the stiffness coefficient *K*eff changes to

$$K\_{\rm eff} = \frac{J\_{\rm piezo}^{\*} \int\_{0}^{L\_{\rm E}} \left(\frac{d^{2}\phi\_{1}(\mathbf{x})}{dx^{2}}\right)^{2} \mathbf{dx} + J^{\*} \int\_{L\_{\rm E}}^{L} \left(\frac{d^{2}\phi\_{1}(\mathbf{x})}{dx^{2}}\right)^{2} \mathbf{dx}}{\phi\_{1}^{2}(\mathbf{x} = L)},\tag{47}$$

for which the term *J*∗ piezo is defined as

$$J\_{\rm piezo}^{\*} = Y\_{\rm s}B\left[\frac{1}{3}\left(\frac{h\_{\rm s} - h\_{\rm P}}{2}\right)^{3} + \frac{1}{3}\left(\frac{h\_{\rm s} + h\_{\rm P}}{2}\right)^{3} - z\_{\rm N}^{2}h\_{\rm s}\right] \\ + Y\_{\rm p}B\left[\frac{1}{3}\left(\frac{h\_{\rm s} + h\_{\rm P}}{2}\right)^{3} - \frac{1}{3}\left(\frac{h\_{\rm s} - h\_{\rm P}}{2}\right)^{3} - z\_{\rm N}^{2}h\_{\rm P}\right] \\ + \nu \,, \tag{48}$$

where

$$\upsilon = \frac{c\_{31}^2}{\epsilon\_{33}^S} B \left[ \frac{1}{3} \left( \frac{h\_{\rm s} + h\_{\rm P}}{2} \right)^3 - \frac{1}{3} \left( \frac{h\_{\rm s} - h\_{\rm P}}{2} \right)^3 - z\_{\rm Pp}^2 h\_{\rm P} \right] \tag{49}$$

and *J* \* is defined as

$$f^\* = \mathcal{Y}\_s B \left[ \frac{1}{3} \left( \frac{h\_\mathsf{s} - h\_\mathsf{P}}{2} \right)^3 + \frac{1}{3} \left( \frac{h\_\mathsf{s} + h\_\mathsf{P}}{2} \right)^3 - z\_\mathsf{N}^2 h\_\mathsf{s} \right] \\ + \mathcal{Y}\_\mathsf{P} B \left[ \frac{1}{3} \left( \frac{h\_\mathsf{s} + h\_\mathsf{P}}{2} \right)^3 - \frac{1}{3} \left( \frac{h\_\mathsf{s} - h\_\mathsf{P}}{2} \right)^3 - z\_\mathsf{N}^2 h\_\mathsf{P} \right]. \tag{50}$$

Next, the damping coefficient *B*eff is defined as

$$B\_{\rm eff} = \frac{2b\_{\rm r}}{\Omega\_1} K\_{\rm eff\ \prime} \tag{51}$$

and the electromechanical coupling coefficient *θ*eff changes to

$$\theta\_{\rm eff} = \frac{\kappa}{\phi\_1(\mathbf{x} = L)} \left. \frac{d\phi\_1}{d\mathbf{x}} \right|\_{\mathbf{x} = L\_{\rm E}} \text{ \,\tag{52}$$

where *κ* is re-calculated with respect to the neutral axis *z*<sup>N</sup> as

$$
\kappa = -\varepsilon\_{31} \left( z\_{\rm TP} - z\_{\rm N} \right). \tag{53}
$$

Then, the effective load *F*eff is defined as

$$F\_{\rm eff} = -\frac{m^\* a\_0 \int\_0^{L-L\_{\rm Mt}} \phi\_1(\mathbf{x}) \mathbf{dx} + \left(m^\* + \frac{M\_0}{L\_{\rm Mt}}\right) a\_0 \int\_{L-L\_{\rm Mt}}^L \phi\_1(\mathbf{x}) \mathbf{dx}}{\phi\_1(\mathbf{x} = L)}\,,\tag{54}$$

and the equivalent capacity *C*eq is defined as

$$\mathcal{L}\_{\text{eq}} = \epsilon\_{33}^{\text{S}} \frac{BL\_{\text{E}}}{h\_{\text{P}}}.\tag{55}$$

#### **3. Verification of Analytical Model Based on Experimental Results**

In this chapter, the derived single DOF model is verified for a time-harmonic kinematic excitation. Three different piezoelectric energy harvesters with known geometry and materials are analyzed and their measured responses are compared with the simulations of the derived single DOF model.

The first experiment is a well-known published work of Erturk and Inman [26] where the authors used a bimorph with PZT-5A piezoelectric material and electrodes spanning over the whole bimorph's length, providing a linear dynamic response. The other two experiments included both bimorph and unimorph configurations of piezoelectric harvesters with a partial electrode length. Both these experiments were conducted in laboratories of Brno University of Technology. The first of these experiments used a bimorph made of PZZN-PLZT piezoceramic [6] which exerted a weak non-linear response in the frequency domain. The second experiment used a simple unimorph configuration for wearables with a thin PVDF layer. Geometrical parameters and material data of individual harvesters for both the piezoelectric layer and the substrate are summarized in Tables 1 and 2, respectively. Data for the PZT-5A bimorph is extracted from [26], the PZZN-PLZT from [6] and the PVDF from [16].


**Table 1.** Parameters of individual harvesters used in experiments.

**Table 2.** Material properties of piezoelectric layers and substrates for each harvester used in experiments.


#### *3.1. PZT-5A Bimorph with a Full Electrode Length and a Linear Response*

This experiment was described and published in detail in paper [26]. This experimental work has a very high impact and for this reason it was used in our analysis as an etalon for the other piezoelectric harvesters. The geometric model of this piezoelectric harvester is in accordance with the model in Figure 1. The bimorph's piezoelectric layers were made of PZT-5A and the substrate was made of brass. It included electrodes covering the whole bimorph's length for harvesting the generated charge. Geometric parameters of the bimorph and properties of used materials are summarized above in Tables 1 and 2, respectively. This bimorph had an experimentally determined damping ratio *b*<sup>r</sup> = 0.027. Since the authors did not state a full description of the tip mass' position and dimensions, it is assumed that the tip mass is located exactly at the bimorph's free end with dimensions allowing for considering the tip mass as a point particle.

This harvester was subjected to a time-harmonic kinematic excitation with a varying forcing frequency *f*. The experiment mapped how amplitudes of generated electrical power and amplitudes of velocity of the bimorph's free end change with a varying forcing frequency upon different values of connected resistive load. Furthermore, the experiment mapped how the peak values of generated electrical power vary with connected resistive load at a specific forcing frequency.

A comparison of published and measured results with the output of our analytical model is presented in Figure 3. The experiment tracked how the peak values of generated electrical power and the velocity amplitude at the beam's free end d*q*0/d*t* change with a varying forcing frequency. The results are displayed for three different values of used resistive load: 1 kΩ, 33 kΩ and 470 kΩ. The graphs show a good match between the output of the analytical model and the obtained experimental data for all three used resistive loads. Note that some discrepancies exist upon the first resonant frequency of the bimorph; this is mainly due to a steep gradient of calculated results near the first resonant frequency. The reader should also note that for resistive loads of 1 kΩ and 33 kΩ there is a slight difference in resonant frequencies between the real bimorph and the analytical single DOF model. This deviation is caused by the used approximative function *φ*<sup>1</sup> which does not account for a concentrated tip mass at the beam's free end. Therefore, as mentioned earlier, the used approximative function forces the beam to behave slightly stiffer and lowers the amount of generated electrical power.

**Figure 3.** Comparison of electrical power and velocity of tip mass for both experimental results [26] and analytical model.

The output of analytical model matches perfectly with the experimental results for both the short-circuit frequency *f* SC and the open-circuit frequency *f* OC of this coupled electromechanical system. The short-circuit and open-circuit frequency are the first resonant frequencies in case of *R*<sup>l</sup> = 0 and *R*<sup>l</sup> → ∞, respectively. The match of simulation results with the measured ones for various values of resistive load *R*<sup>l</sup> and kinematic excitation at both the short-circuit frequency *f* SC and the open-circuit frequency *f* OC is shown in Figure 4. Both states correspond with operations slightly below and above the resonance excitation for various values of resistive load, which determine the value of actual resonance frequency.

**Figure 4.** Peak power values as a function of resistive load upon excitation at short-circuit resonance frequency and the open-circuit resonance frequency.

While in case of the open-circuit forcing frequency the results of the analytical single DOF model agree with the measured values, for the short-circuit case the calculated values from the single DOF model are slightly shifted towards higher values of resistive load. Also note that there are differences in both frequencies between the real bimorph and the single DOF model. While in case of the open-circuit frequency this difference is very small, for the short-circuit frequency this difference is notably larger and affects the value of optimal resistive load for which the generated electrical power reaches its peak value. This is caused by the used approximative function *φ*<sup>1</sup> which causes the beam to behave stiffer. Nevertheless, this inaccuracy is negligible in terms of using analytical models for a rough prediction of the generated power when the system is excited by real vibrations.

#### *3.2. PZNN-PLZT Bimorph with Partial Electrode Length and Weak Non-Linear Response*

The experiment with PZNN-PLZT bimorph (see Figure 5) was conducted in a laboratory at Brno University of Technology with a partial electrode length (there are no electrodes under the tip mass). Geometrical parameters of this bimorph are summarized above in Table 1. The bimorph's piezoelectric layers were made of PZNN-PLZT, which is in detail described in [6], and the substrate was made of a common steel shim. Properties of these materials are listed above in Table 2. Electrodes were made using a thin silver tape casting. Since the silver electrodes were substantially thinner than other layers, they were not accounted for in the calculation of single DOF model parameters due to their negligible effect on the net mass and the beam's stiffness. The bimorph had an experimentally determined damping ratio *b*r = 0.025 via an impulse response in the short-circuit state.

The clamping of the used bimorph was kinematically excited at several forcing frequencies near the bimorph's first resonant frequency with a constant acceleration amplitude *a*<sup>0</sup> = 0.1 g. The aim of this experiment was to track results, namely the RMS of generated voltage and RMS of velocity of the tip mass, at different excitation frequencies close to the bimorph's first natural frequency. Also, the optimal resistive load was sought at which the bimorph generates maximal electrical power at its current first resonant frequency which slightly varies with changes in *R*l.

A comparison between the measured data and the calculated output of the analytical model is shown in Figure 6, namely the RMS values of output voltage and those of velocity of the bimorph's free end as a function of forcing frequency. The results are displayed for two values of used resistive load: 1 MΩ and 10 MΩ. The measured data shows a weak non-linear softening dynamic behavior; however, the analytical single DOF model with linearized parameters still shows a very good degree of accuracy for both used resistive loads in terms of achieved amplitudes.

**Figure 5.** PZNN-PLZT bimorph used in experiment.

**Figure 6.** Comparison of generated voltage and velocity of harvester's tip mass obtained from measurement and developed analytical model for *R*<sup>l</sup> are 1 MΩ and 10 MΩ.

Our experiment also tracked the values of generated electrical power as a function of used resistive load. During the measurement, the forcing frequency was adjusted for each value of resistive load so that it matched the bimorph's actual first resonant frequency. Both the analytical model and the experiment show (Figure 7) that the optimal resistive load is approx. 1.5 MΩ and, at the same time, also the maximal values of generated electrical power calculated with the single DOF model agree with experimental data at

all used values of resistive load. The reader should note here that the curve from the analytical model is slightly shifted to higher values of *R*<sup>l</sup> which is, similarly as in the previous experiment, due to the used approximative function *φ*1, which makes the beam model behave slightly stiffer.

**Figure 7.** Power as function of resistive load upon excitation at actual resonant frequencies.

#### *3.3. PVDF Unimorph with a Partial Electrode Length and a Linear Response*

PVDF piezoelectric energy harvesters are very often presented as a suitable kinetic energy harvester [29] and for this reason the PVDF material was chosen for the last experiment, which was also conducted in a laboratory at Brno University of Technology. The PVDF foil is used in a unimorph configuration of a clamped cantilever with a partial electrode length shown in Figure 8. Parameters of this unimorph are listed above in Table 1. The unimorph's piezoelectric layer is a PVDF foil and the substrate is a steel shim [30]. Properties of these materials are summarized above in Table 2. Electrodes were made using a thin silver tape casting. The silver electrodes were not accounted in the calculation of single DOF model parameters as in the previous model due to their negligible effect on the net mass and beam's stiffness. The clamping of the used unimorph was kinematically excited at several forcing frequencies near the unimorph's first natural frequency (*f* 1,r = 18.7 Hz) with a constant acceleration amplitude *a*<sup>0</sup> = 0.035 g. The unimorph had an experimentally determined damping ratio *b*r = 0.0065 via an analysis of impulse response in the short-circuit state.

This experiment measured the RMS of output voltage and the amplitude of velocity of the tip mass at different forcing frequencies close to the unimorph's first natural frequency. A comparison between the measured data and the calculated output of the analytical model is shown in Figure 9, namely the RMS values of output voltage *U* and amplitudes of velocity of the tip mass as a function of a forcing frequency for *R*<sup>l</sup> = 10 MΩ. One can see that the first resonant frequency of the analytical model is again slightly higher due to used approximative function *φ*1; nevertheless, the calculated values from the analytical model agree with the measured ones.

#### *3.4. Single DOF Model Parameters of Considered Harvesters*

The calculated parameters of individual harvesters which were used as input in analytical models are summarized in Table 3. The values of effective load *F*eff were normalized with respect to 1 g of base acceleration.

**Figure 8.** PVDF unimorph used in the experiment and a strip of PVDF foil used as the piezoelectric layer.

**Figure 9.** Comparison of output voltage and velocity of tip mass obtained from the measurement and by using the developed analytical model for *R*<sup>l</sup> = 10 MΩ.

**Table 3.** Parameters of each piezoelectric harvester used in the analytical model.


The steady-state results calculated with the developed analytical single DOF model and parameters given in Table 3 showed an excellent agreement with all presented experiments. Combined with low usage of computer resources, the developed analytical model presents a simple and very effective tool for proper designing of piezoelectric harvesters. Moreover, the model can be used to simulate transient responses of the considered harvester (represented by its geometry and materials) to arbitrary time-dependent loads as demonstrated further in the text. For increased accuracy outside the vicinity of the

harvester's first resonant frequency, additional mode shapes (*φ*2, *φ*3, etc.) can be supplemented. Moreover, the single DOF model can easily be extended to support calculations of strain and stress levels within the beam's layers to determine a maximal allowable load as shown in [31].

#### **4. Comparison of Piezoelectric Materials for Kinetic Energy Harvesting Purposes**

Since the analytical single DOF model of a piezoelectric harvester was successfully validated for both unimorph and bimorph configurations using the data from three different piezoelectric materials and experiments, the verified models of three harvesters will enable to determine the effectivity of used piezoelectric materials in different energy harvesting applications. Here, the three materials considered in the scope of this work (PZT-5A, PZZN-PLZT and PVDF) are compared in terms of harvested electrical power when subjected to harmonic vibrations (lab shaker) and in terms of harvested electrical energy when subjected to random vibrations (human body movement).

#### *4.1. Harmonic Vibrations Case*

To compare the output of harmonically excited piezoelectric harvesters made of different piezoelectric materials, their dynamic parameters must be similar, that is, their effective mass *M*eff and eigenfrequency *f*1. This can be done by changing dimensions of the considered harvesters; however, doing this will also lead to changes in the piezoelectric coupling coefficient and ultimately making the comparison invalid. To overcome this issue and maintain comparability, the volume of polarized piezoelectric materials was kept constant. The dimensions of polarized piezoelectric material in case of piezoceramic bimorphs were fixed at values *L*<sup>E</sup> × *B* × *h*<sup>p</sup> = 40 × 10 × 0.26 mm and in case of PVDF unimorph at *L*<sup>E</sup> × *B* × *h*<sup>p</sup> = 40 × 40 × 0.13 mm due to manufacturing limits of PVDF foils (these must be thin but can span over a large area). Then, to reduce the complexity of this optimizing task, the value of *L*Mt was fixed at 5 mm and the thickness of the substrate *h*<sup>s</sup> was fixed at 0.15 mm (0.3 mm) in case of bimorphs (PVDF unimorph). Thus, the only parameters left for optimizing were the total length of the harvester *L* and and the tip mass *M*t. These parameters were then tuned (see Table 4) to achieve values of *M*eff and *f*<sup>1</sup> common to all three harvesters. Upon the study, the harvesters were forced with a time-harmonic base acceleration for various values of resistive load *R*<sup>l</sup> at their actual resonant frequencies.


**Table 4.** Tuned dimensions of harvesters and their equivalent single DOF model parameters used in the comparison.

The comparison (see Figure 10) revealed that PZT-5A is the best suitable material of the considered ones for energy harvesting purposes thanks to its high power output (several mW per 1 g of base acceleration) and a broad range of optimal resistive load (approx. 150 kΩ to 1.1 MΩ) due to its strong piezoelectric coupling. PZZN-PLZT is also suitable for energy harvesting applications since it offers high power output of about 1 mW per 1 g of base acceleration for resistive loads close to 1 MΩ. On the other hand, PVDF generates the least amount of power of the three materials and due to its very high optimal resistive load it is not sufficient for energy harvesting applications. Note that the strong piezoelectric coupling in case of PZT-5A harvester significantly damps the power output between the two optimal resistive loads which results in a local minimum surrounded by two local maxima.

**Figure 10.** Comparison of harvested electrical power among the considered materials for various resistive loads upon simple harmonic forcing at the first resonant frequency.

#### *4.2. Random Vibrations Case*

For this comparison, typical mechanical vibrations of a human forearm which are encountered in wearables applications were measured [32] and analyzed using the developed model. The measured time-course of acceleration *a*(*t*), see Figure 11a, is generated by a random movement of the forearm and can be thought of as a representative of random vibrations as can be seen from its spectrogram in Figure 11b. Therefore, it is perfect for the comparison of energy harvesting devices since a steady-state response, whose magnitude depends on how close the forcing frequency is to the harvester's resonant frequency, will not occur.

The measured acceleration *a*(*t*) is used as input for a transient analysis of the derived single DOF analytical model, where the applied force is a function of acceleration data. This model of coupled electro-mechanical system is realized in Matlab Simulink simulation environment [33] and its aim is to track the amount of harvested electrical energy for various values of resistive load *R*l. Simulation results of predicted harvested energy for this wearable operation for the harvesters used in the experiments and the tuned harvesters from Section 4.1 are shown in Figure 11c,d, respectively.

In case of unmodified harvesters used in the experiments (Figure 11c), the PZT-5A harvester is able to convert the most of mechanical energy (~0.17 mJ) among the three compared harvesters and has a low value of optimal resistive load (~118 kΩ). The energy output of PZZN-PLZT harvester (~0.07 mJ) is of the same order as the one of the PZT-5A harvester, however its much higher optimal resistive load (~1.82 MΩ) makes it less suitable for energy harvesting applications. On the contrary, the PVDF harvester is not suitable for energy harvesting applications at all due to very low energy output (~0.07 μJ), but it will find its use in sensing applications due to very high value of optimal resistive load (~28 MΩ). The same also applies to results in case of tuned harvesters (Figure 11d), where the tuned PZT-5A harvester once again shows that the energy harvesting properties of PZT-5A are far superior to those of PZZN-PLZT and PVDF. The increase in harvested

electrical energy for the tuned PZT-5A harvester compared to the unmodified geometry is due to its lower resonant frequency which was reduced from 46.8 Hz to 31.1 Hz.

**Figure 11.** (**a**) Measured acceleration of a random movement of a human wearable and (**b**) its spectrogram; a comparison of harvested electrical energy from the human forearm movement among (**c**) the piezoelectric harvesters used in the experiments and (**d**) the tuned piezoelectric harvesters from Section 4.1 for different values of resistive load.

The results of both these comparisons showed that both PZT-5A and PZZN-PLZT piezoceramic harvesters are suitable for energy harvesting purposes, although operating at different values of resistive load, whereas the harvester based on piezoelectric polymer PVDF provides an insufficient energy harvesting system due to a very low amount of harvested energy. However this energy harvester should primarily be used in sensing applications [34].

#### **5. Conclusions**

The main aim of this paper was to compare the effectiveness of materials commonly used in energy harvesting operations using a single DOF model and at the same time analyze the effect of used mode shape function on simulation results. The single DOF model of a cantilever piezoelectric harvester in both bimorph and unimorph configurations was derived based on Euler-Bernoulli beam theory. Output of the model was confronted with available experimental data obtained from three different piezoelectric harvesters (PZT-5A bimorph, PZZN-PLZT bimorph and PVDF unimorph) and showed a good degree of accuracy. It is obvious that the presented model of an energy harvester can be used for various piezoelectric materials. Therefore, the developed single DOF analytical model represents a simple and very helpful tool for designing piezoceramic vibration energy harvesters. Moreover, it could easily be employed to check if a particular kinetic energy harvester provides sufficient output power for the intended application. Or inversely, the model could be used to design a piezoceramic harvester with optimized operational parameters and dimensions due to the model's ability to predict the amount of harvested energy in particular operational conditions. Moreover, the developed model can easily be extended to support calculations of strain and stress levels within harvester's layers for further assessments concerning strength and fatigue limits.

The model is primarily intended for operations of the harvester at frequencies where vibrations consist mostly of the first mode shape, since it offers the best operational conditions for energy harvesting (no strain nodes). If higher vibrational modes are of interest, the developed model can easily be extended by supplementing their respective shape functions and using the superposition principle. It was found that the quality of the used approximative function for the first mode shape affects the beam model's stiffness in such a way so that simpler (less accurate) approximative functions force the beam model to behave stiffer, i.e., its resonant frequency being shifted to higher values. Also, the way how device layers and electrodes are assembled can also affect the stiffness of the system and could potentially result in a weak nonlinearity, which was observed in one of our experiments. Nevertheless, a typical assembly of layers and the approximative shape function used in this work (the shape of the beam's first vibrational mode without a tip mass) still shows a good degree of accuracy.

The single DOF model itself poses as a very effective tool whose main advantages are low computer resources usage and the ability to calculate transient responses for arbitrary time-dependent loads. Both these features were employed in an energy harvesting effectivity comparison of the three materials used in the scope of this work. The materials comprised of PZT-5A, which is known nowadays to be one of the best materials for energy harvesting purposes, and PZZN-PLZT and PVDF which are used in our laboratory for designing vibration energy harvesting devices. The comparison was split into two parts with respect to forcing: case of simple harmonic vibrations and case of random vibrations. The case of simple harmonic vibrations was carried out so that the harvesters' dimensions were tuned in order to achieve common value of seismic mass and resonant frequency for all three harvesters and at the same time the harvesters had same volume of polarized piezoelectric material. In case of random vibrations the harvesters were subjected to a non-harmonic and non-periodic vibrations typical for wearables applications. Results from both comparison cases showed that the piezoceramic harvesters (PZT-5A and PZZN-PLZT) are a perfect choice for energy harvesting applications, though geometry and electrical load must be optimized. On the contrary, the PVDF harvester is not suitable for energy

harvesting purposes due to very low values of harvested energy despite many recent papers reporting the otherwise, and its potential lies in sensing applications.

**Author Contributions:** Conceptualization, Z.H.; methodology, Z.M. and O.S.; software, Z.M. and O.S.; validation, O.R.; formal analysis, Z.M. and O.S.; investigation, Z.M. and O.R.; resources, Z.H.; data curation, O.R.; writing—original draft preparation, Z.M.; writing—review and editing, Z.H. and O.S.; visualization, Z.M.; supervision, Z.H. and O.S.; project administration, Z.H.; funding acquisition, Z.H. All authors have read and agreed to the published version of the manuscript.

**Funding:** This presented research and development was supported by the Czech Science Foundation project GA19-17457S "Manufacturing and analysis of flexible piezoelectric layers for smart engineerin", Czech Republic.

**Institutional Review Board Statement:** Not applicable.

**Informed Consent Statement:** Not applicable.

**Conflicts of Interest:** The authors declare no conflict of interest. The funders had no role in the design of the study; in the collection, analyses, or interpretation of data; in the writing of the manuscript, or in the decision to publish the results.

#### **References**


## *Article* **A New Method to Perform Direct Efficiency Measurement and Power Flow Analysis in Vibration Energy Harvesters**

**Jan Kunz 1, Jiri Fialka 1,2, Stanislav Pikula 1, Petr Benes 1,2,\*, Jakub Krejci 1, Stanislav Klusacek 1,2 and Zdenek Havranek 1,2**


**Abstract:** Measuring the efficiency of piezo energy harvesters (PEHs) according to the definition constitutes a challenging task. The power consumption is often established in a simplified manner, by ignoring the mechanical losses and focusing exclusively on the mechanical power of the PEH. Generally, the input power is calculated from the PEH's parameters. To improve the procedure, we have designed a method exploiting a measurement system that can directly establish the definitionbased efficiency for different vibration amplitudes, frequencies, and resistance loads. Importantly, the parameters of the PEH need not be known. The input power is determined from the vibration source; therefore, the method is suitable for comparing different types of PEHs. The novel system exhibits a combined absolute uncertainty of less than 0.5% and allows quantifying the losses. The approach was tested with two commercially available PEHs, namely, a lead zirconate titanate (PZT) MIDE PPA-1011 and a polyvinylidene fluoride (PVDF) TE LDTM-028K. To facilitate comparison with the proposed efficiency, we calculated and measured the quantity also by using one of the standard options (simplified efficiency). The standard concept yields higher values, especially in PVDFs. The difference arises from the device's low stiffness, which produces high displacement that is proportional to the losses. Simultaneously, the insufficient stiffness markedly reduces the PEH's mechanical power. This effect cannot be detected via the standard techniques. We identified the main sources of loss in the damping of the movement by the surrounding air and thermal losses. The latter source is caused by internal and interlayer friction.

**Keywords:** piezoelectric; piezoelectric ceramic; lead zirconate titanate (PZT); polyvinylidene fluoride (PVDF); energy harvesting; efficiency; efficiency measurement; power conversion; power flow

#### **1. Introduction**

The increasing demands on the monitoring and predictive maintenance of mechanical structures, intelligent building control, and other remote sensing applications generate the need to install electronics at isolated locations. Such sites are often characterized by difficult access and a lack of infrastructure. Self-powered electronic devices thus became an essential prerequisite. These instruments and apparatuses gather electricity from sources such as light, temperature difference, and diverse forms of kinetic energy, including ambient mechanical vibrations [1,2]. The vibrations can be converted to usable electrical energy by means of piezoelectric energy harvesters (PEHs).

This type of movement is present in machines and structures where collecting relevant operational information has a beneficial impact, namely in industrial equipment, cars, aircraft, and bridges [1,3]. Moreover, PEHs are also employed in harvesting from human

**Citation:** Kunz, J.; Fialka, J.; Pikula, S.; Benes, P.; Krejci, J.; Klusacek, S.; Havranek, Z. A New Method to Perform Direct Efficiency Measurement and Power Flow Analysis in Vibration Energy Harvesters. *Sensors* **2021**, *21*, 2388. https://doi.org/10.3390/s21072388

Academic Editor: Amir H. Alavi

Received: 5 March 2021 Accepted: 25 March 2021 Published: 30 March 2021

**Publisher's Note:** MDPI stays neutral with regard to jurisdictional claims in published maps and institutional affiliations.

**Copyright:** © 2021 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (https:// creativecommons.org/licenses/by/ 4.0/).

movement [1]; substantial effort has therefore been invested in developing and improving PEHs to generate as much energy as possible from the oscillations.

The power output constitutes the key criterion to determine the applicability of a PEH; however, this factor does not provide sufficient information if individual PEHs are to be compared. The reason is that the maximum power output is defined mainly by the PEH's dimensions and vibration amplitude. In view of these facts, Beeby et al. [4] proposed normalized power density (NPD), which divides the power output by the vibration amplitude and the volume of the harvester. Although such metrics appear suitable for comparing PEHs, the authors also mention certain drawbacks. One of the disadvantages rests in ignoring important factors, including the bandwidth [4]. Generally, however, Beeby et al.'s approach is simple, straightforward, and easily applicable. Other researchers developed more sophisticated methodologies, involving the frequency parameters; these techniques were analyzed comprehensively by Hadas et al. [5], who outlined their benefits and drawbacks.

A portion of relevant papers interpret the PEH as an energy converter, meaning that the power output and related alternatives are not considered convenient for a comparison of PEHs. The authors thus employ efficiency as the criterion to examine the devices, as is common in other energy converters. The relevant theories and measurement procedures were summarized and discussed by Yang et al. [6]. In their paper [6] the efficiency not only embodies an essential factor in the development and optimization of PEHs but also finds use as a parameter allowing comparison between energy harvesting methods.

#### *Efficiency*

The efficiency is the ratio between the power output and input. In PEHs, the electrical power output is measurable without difficulty. Determining the mechanical power input, however, constitutes a more challenging task, as the process cannot be easily monitored [6]. To establish the efficiency, the power input of the harvester is often simplified to its mechanical power [6,7] because this quantity can be easily calculated via the PEH's parameters.

The calculation is performed by means of, for example, the widely applied formula proposed by Richards et al. [7], which defines the efficiency for the resonant frequency and optimal load [6]:

$$\eta = \frac{1}{2} \frac{k\_{\rm sys}^2}{1 - k\_{\rm sys}^2} \bigg/ \left(\frac{1}{Q} + \frac{1}{2} \frac{k\_{\rm sys}^2}{1 - k\_{\rm sys}^2}\right) \tag{1}$$

where *η* is the efficiency, and *k*sys and *Q* denote the system coupling coefficient and the quality factor of the PEH, respectively.

Some of the researchers [6,8,9] expanded the efficiency calculation to involve broader frequency and load ranges. Regrettably, neither the calculating process nor the measurement of the parameters are unified, resulting in that the determined efficiencies oscillate between less than 1% [10,11] and more than 80% [8,12] in lead zirconate titanate (PZT) PEHs.

The role of the system coupling coefficient (*k*sys) is not identical with that of the electromechanical coupling coefficient (*k*eff) defined in [13,14]. The latter embodies a piezomaterial coefficient, whereas the former (*k*sys) describes the whole structure of the PEH. The material coupling coefficient can be obtained from the material resonance and anti-resonance frequencies [13]:

$$k\_{\rm sys}^2 = \frac{f\_{\rm o}^2 - f\_{\rm s}^2}{f\_{\rm o}^2},\tag{2}$$

where *f*o and *f*s are the open and the short natural frequencies of the PEH, respectively.

Equation (2) has found wide use [6–8], in spite of the fact that it is valid only for undamped systems (damping coefficient, *ζ* = 0) [15,16]. Utilizing Formula (2) in damped systems (*ζ* > 0) produces inaccurate parameters [15] and, consequently, imprecise efficiency.

Regrettably, this method (Equation (2)) for calculating the *k*sys exhibits significant uncertainty due to two large numbers being subtracted close to each other; therefore, the efficiency calculation (Equation (1)) will comprise a major error, too. This error, however, is not evaluated in the literature [6,7,10,11,17,18].

Some experts [6,12] also measured experimentally the mechanical power of PEHs. Although such an attitude finds application in validating the derived formulas, it still ignores the mechanical losses. By extension, the parameters of the devices have to be employed to calculate the power. This simplification assumes that the ambient vibration source has an "infinite power" compared to the consumed power of the PEH. The mechanical losses can therefore be ignored, as they are easily replenished from the source [19]. In some situations, such as harvesting from human motion, this assumption is invalid, and the mechanical losses have to be taken into account [19]. Furthermore, the mechanical losses also constitute an important parameter in developing the PEH.

Liao and Sodano [20] established that the approach neglecting the mechanical losses is oversimplified, and they proposed that the input power should not be considered the mechanical power of the PEH. The authors recommend using the power fed into the device from the ambient source to maintain the steady state vibrations. Although this claim is correct, the calculation of the power from the parameters of the PEH cannot be characterized in the same manner. The reason rests in that the procedure involves only the internal and electrical losses and ignores the mechanical ones (e.g., the internal friction or aerodynamic drag). In this context, Yang et al. [6], for example, found close similarity between their method and the PEH's damping coefficient.

Based on our review of the literature, we can conclude that the mechanical losses are disregarded in the current efficiency measurements, which exploit the mechanical power of PEHs as the input. Furthermore, although the harvesters' parameters are utilized for the calculations, the measurement method and the uncertainties are either characterized unclearly or remain ignored [6,8,10,11].

Such issues then justify the need for a new measurement system capable of evaluating the efficiency via the mechanical power supplied to a PEH by an ambient vibration source (proposed efficiency). This power is determined from the force and velocity of the source (in our case, a vibration shaker), eliminating the necessity to use the parameters of the PEH in the calculation. Furthermore, the same measurement assesses the harvester's mechanical power, too, and the efficiency can then be established in a simplified manner (simplified measured). Importantly, the device's parameters (*k*sys and *Q*) are calculable from the results to facilitate computing the efficiency according to (1) (simplified calculated). In addition, the uncertainties of all of the aforementioned approaches are determined during the measuring procedure.

As indicated, the unique setup allows us to yield the proposed efficiency, the simplified efficiencies, and the related uncertainties from merely a single measurement. The methods are thus compared in identical conditions. Expectably, the simplified efficiencies, being comparable, should exhibit the same values. The proposed efficiency nevertheless shows lower values, as it includes the mechanical losses.

The paper is organized as follows: The introductory sections set out the problem, together with relevant processes, procedures, and options. Section 2 debates the methodology, characterizing the power flow in piezoelectric harvesters, outlining the calculation of the input power, and describing the measurement system and its uncertainties. In Section 3, the main discussion and results are presented, including the outcomes obtained from a measurement of the efficiency of two different commercially available harvesters. The results of the proposed method are correlated to common simplified efficiencies. Finally, the power flow is discussed, in conjunction with a quantitative analysis of the losses.

#### **2. Methods**

Prior to describing the efficiency measurement system, we will discuss the power flow in PEHs (Figure 1). During each cycle, the energy is transmitted back and forth due to the oscillations of the system. In this chapter, however, only the power flow average through several cycles is assumed.

**Figure 1.** The power flow in a piezo energy harvester (PEH). The traditional approach to determine the efficiency considers the harvester's mechanical power (*P*M) as the input power. Alternatively, according to the definition, the power fed into the PEH from the ambient source (*P*VIB) has to be used as the input.

In a laboratory environment, the vibrations are generated by a vibration exciter (shaker), which produces a vibration power (*P*SHPEH ) and has its own losses (*P*SHNO ).

The power from the vibration source (*P*VIB) is transferred to the vibration movement of the PEH, i.e., the harvester's mechanical power (*P*M). Similarly to every other power transfer, this process generates some losses. In the presented case, the two main sources are the thermal losses (*P*TH) generated by the internal friction in the PEH, and the acoustic losses (*P*AIR), originating from the movement of the PEH in the air.

The mechanical power (*P*M) is further transferred to the harvester's usable electrical power output (*P*E). Here, the main sources of loss are the internal losses (*P*IL), arising mainly from the electromechanical losses in the piezoelectric and electrical dissipation (*P*EL) on the PEH's internal resistance.

In steady state vibrations, however, only a part of the PEH's mechanical power (Δ*P*M) is dissipated in the form of the internal and electrical losses. The electrical power output is replenished from the vibration source, as proposed by [20]. Nevertheless, the amount of such replenished power is difficult to evaluate. For simplification we assume in the calculation that the entire mechanical power of the PEH is replenished from the vibration source. This approach also allows us to be consistent with other efficiency-related papers, except for [20]. Similarly to our procedure, however, the authors of this referenced source, too, were unable to measure the losses.

The applied simplification underestimates the thermal losses (*P*TH) as these are calculated via the PEH's mechanical power (*P*M). The efficiencies are not affected, because the traditional approach to determine the efficiency employs the mechanical power (*P*M) as the input power [6]. However, according to the definition [20], the power supplied to the PEH from the ambient source, namely the *P*VIB in our case, has to be used as the input. For clarification purposes, the methods to measure or calculate the powers expressed in Figure 1 are described in the sections below.

#### *2.1. Shaker Vibration Power*

The vibration power of the shaker (*P*SH) is calculated out of the known force and velocity (Equation (3)) [6,18,21]. As both of the quantities are vectors, their phases have to be measured, too. In the steady state condition and at harmonic vibration, the power of the mechanical vibration is calculated as:

$$p\_{\rm SHI} = \mathbf{F} \cdot \boldsymbol{\mathfrak{v}} = |\mathbf{F}| |\boldsymbol{\mathfrak{v}}| \cos(\boldsymbol{\phi}).\tag{3}$$

The velocity of the structure's base is either directly measured or, as in our case, calculated from the known acceleration. Since the vibrations are steady state and harmonic, the calculation simplifies into:

$$|\mathbf{v}| = \left| \int \mathbf{a} dt \right| = \int a\_{\rm ampl} \sin(2\pi f t + \phi\_a) dt = \frac{a\_{\rm ampl}}{2\pi f} \cos(2\pi f t + \phi\_a). \tag{4}$$

The electrodynamic shaker is, principally, a loudspeaker. To calculate the force generated by the shaker, we can then employ the formula of Ampere's law, commonly used for this purpose in acoustics [22]:

$$F = lI \times \mathbf{B}.\tag{5}$$

As the magnetic flux density *B* in the coil gap and the length of the wire *l* in the magnetic field remain constant, their product, *kBl*, can be easily determined [22], simplifying Equation (5) to read

$$|F| = k\_{Bl}|I| = k\_{Bl}I\_{\text{ampl}}\sin(2\pi ft + \phi\_l). \tag{6}$$

Substituting (4) and (6) into (3) and averaging the outcome through the period, we yield:

$$\begin{split} P\_{\rm SH} = \frac{1}{T} \int\_{T} p\_{\rm SH} dt = \frac{1}{T} \int\_{T} k\_{\rm BI} I\_{\rm amp} \sin(2\pi f t + \phi\_{l}) \frac{a\_{\rm ampl}}{2\pi f} \cos(2\pi f t + \phi\_{a}) dt \\ &= k\_{\rm BI} I\_{\rm rms} \frac{a\_{\rm rms}}{2\pi f} \cos\left(\phi - \frac{\pi}{2}\right), \quad (7) \end{split} \tag{7}$$

where *p*SH and *P*SH are the shaker's instantaneous and average vibration power, respectively; *F*rms and *v*rms are the effective values of the vibration force and velocity, respectively; *kBl* denotes the shaker's Bl coefficient; *I*rms and *a*rms are the effective values of the drive current and acceleration, respectively; *φ* is the phase between those two vectors; and *f* represents the frequency of the vibrations.

#### *2.2. Vibration Power*

The measurement of the vibration power (*P*VIB) delivered from the vibration source to the PEH is rather complicated [6]. In the case of steady state harmonic vibrations, however, the power can be determined from the vibration power of the shaker (*P*SH).

The power generated by the vibration shaker (*P*SH) depends on the vibration amplitude, frequency, and load on the vibration rod and fixture. Advantageously, the power to vibrate structures other than the PEH is measurable as well, by the same measurement but without the PEH mounted. Then, the vibration power (*P*VIB) can be calculated simply by subtracting the vibration power without the PEH (*P*SHNO ) from that with the device (*P*SHPEH ). We have:

$$P\_{\rm VII} = P\_{\rm SFIPE} - P\_{\rm SFINO'} \tag{8}$$

where *P*VIB is the vibration power consumed by the PEH from the ambient vibration source, and *P*SHPEH and *P*SHNO stand for the shaker's vibration power with and without the PEH mounted, respectively.

#### *2.3. Acoustic Losses*

The motion of the piezoelectric cantilever is damped by the ambient air. The damping produces an acoustic power, calculable from the known deflection of the moving object. The acoustic power for a harmonic oscillation is given by [22], as follows:

$$P\_{\rm air} = AI = A \frac{1}{2} \rho c \omega^2 d^2 \,, \tag{9}$$

where *A* is the area, *I* denotes the sound intensity, *ρ* represents the air density, *c* is the speed of sound in the air, *ω* stands for the angular frequency, and *d* is the oscillation amplitude.

Due to the nature of the cantilever deflection, the calculation can be materialized when the resonance mode is measured at *n* points evenly distributed on the surface of the cantilever. The measured deflections (*di*) are applicable in establishing the partial intensities (*Ii*), and we sum them to estimate the total acoustic power; we then have:

$$P\_{air\_{tot}} = \sum\_{i=1}^{n} \frac{A}{n} I\_n = \sum\_{i=1}^{n} \frac{A}{n} \frac{1}{2} \rho c \omega^2 d\_i^2. \tag{10}$$

When the deflections are measured on the edge of the structure, this arrangement has to be taken into account and then reflected in the calculation by an appropriate reduction of the corresponding area.

#### *2.4. Thermal Losses*

The power transfer from the vibration rig to the movement of the PEH, the PEH's mechanical power (*P*M), involves losses, such as the internal friction between and in the material layers of the harvester. The losses dissipate some of the power to heat. By extension, in more general terms, losses also accompany the mechanical–electrical power conversion in the piezoelectric (Section 2.8) and the electrical losses (Section 2.7) on the internal resistance; both are also dissipated as heat. Thus, the powers of the individual loss components cannot be evaluated by measuring the heat.

Moreover, in PEHs, the overall power dissipated into heat is so small that it produces only an almost insignificant change in the harvester's surface temperature, for instance, in the PZT PEH at 1 grms the warming of the PEH caused by the thermal losses (14.93 mW; Table 1) is less than 0.02 K. For this reason, we were unable to measure the outcome.

Beneficially, the amount of thermal losses (*P*TH) can be determined in an indirect manner, via:

$$P\_{\rm TH} = P\_{\rm VIB} - P\_{\rm AIR} - P\_{\rm M} \tag{11}$$

where *P*TH and *P*AIR are the thermal and the acoustic losses, respectively, and *P*VIB and *P*<sup>M</sup> denote the vibration and the mechanical powers of the PEH, respectively.

#### *2.5. PEH Mechanical Power*

The mechanical power of a PEH (*P*M) can be measured via Formula (12), derived by Yang et al. [6]:

$$P\_{\rm M} = \frac{1}{2} m\_{\rm 6} a\_{\rm \chi} v\_{\rm Z} \sin(\phi\_{\rm X}) , \tag{12}$$

where *m*e is the PEH's effective mass, *a*x denotes the amplitude of the relative acceleration of the PEH's tip, *v*<sup>Z</sup> represents the amplitude of the base velocity, and *φ*<sup>X</sup> is the phase difference between the base velocity and the PEH's tip acceleration.

#### *2.6. Electrical Power Output*

The electrical power output of the PEH (*P*E) is measured on the resistive load; however, it alternates with the vibration frequency. Thus, the output has to be rectified to be able to power some types of electronics, and additional losses occur. For this reason, rectification

was not involved in our experiments. Relevant power management circuits are also being intensively developed [5]; these are nevertheless not discussed in this paper. We have:

$$P\_{\rm E} = \frac{\mathcal{U}\_{\rm rms}^2}{R\_{\rm L}} \prime \tag{13}$$

where *P*<sup>E</sup> is the PEH's electrical power output, and *U*rms denotes the effective value of the voltage generated by the PEH on the preset resistive load, *R*L.

#### *2.7. Electrical Losses*

To ensure the maximum flow of power from the PEH to the load, the optimum external load has to have the same value as the internal resistance of the PEH. Then, the power transfer is at a maximum but provides an efficiency of only 50%, meaning that half of the generated electrical power is the electrical power output (*P*E) and the other half comprises the electrical losses (*P*EL) dissipated on the internal resistance. Then, the electrical losses (*P*EL) can be easily determined:

$$P\_{\rm EL.} = P\_{\rm E.} \tag{14}$$

where *P*EL denotes the electrical losses on the internal resistance, and *P*<sup>E</sup> is the electrical power output.

#### *2.8. Internal Losses*

The mechanical–electrical power conversion is performed by the piezoelectric layer in the PEH. The portion of the mechanical power that is absorbed by the layer but remains unconverted into the electrical power (*P*E) forms the internal losses (*P*IL), which are dissipated into heat.

The efficiency of the conversion in the piezomaterial is defined by the electromechanical coupling coefficient (*k*eff). The relevant value usually varies between 0.4 and 0.6 in the PZT ceramic but amounts to less than 0.3 in the polyvinylidene fluoride (PVDF) film. The conversion efficiency of the material enclosed in the PEH's structure is slightly different, and therefore the internal losses (*P*IL) cannot be calculated directly by using this coefficient. However, indirect calculation is possible; we then have:

$$P\_{\rm IL} = P\_{\rm M} - P\_{\rm E} - P\_{\rm ELJ} \tag{15}$$

where *P*IL and *P*EL are the internal and electrical losses, respectively, and *P*<sup>M</sup> and *P*<sup>E</sup> represent the mechanical and electrical powers, respectively.

#### *2.9. Efficiency Calculation*

To include the mechanical losses in the efficiency, we consider as an input power the vibration power supplied to the PEH from an ambient vibration source, *P*VIB. We have:

$$\eta = \frac{P\_{\rm E}}{P\_{\rm VIB}} = \frac{\frac{\bar{U}\_{\rm rma}^2}{R\_{\rm L}}}{P\_{\rm SF\_{\rm PEH}} - P\_{\rm SFhNO}},$$

where *η* stands for the measured efficiency; *P*VIB is the PEH's power input; *P*<sup>E</sup> denotes the electrical power output, i.e., the voltage *U* generated across the load resistor, *R*L; and *P*SHPEH and *P*SHNO represent the vibration power with and without the PEH mounted on the shaker, respectively.

#### *2.10. Efficiency Measurement*

To measure the efficiency, we improved the automated measurement system [23] to manage the vibration amplitudes at a lower error rate and to reduce the noise. Moreover, the data post-processing stage was redesigned to compute the power consumption of the PEH. The actual setup is described in detail within Section 2.12.

Because the efficiency results appear to be the most interesting near the PEH's mechanical resonance, this resonance is localized by using sweep sine vibration.

To obtain proper efficiency results, we need to ensure that the PEH vibrates in its first longitudinal mode without any torsional oscillation or clamping. Such a condition was verified via a mode measurement system utilizing a scanning laser vibrometer. This step is especially significant when the examined PEH has unknown parameters. Importantly, the measured mode can also be used for calculating the aerodynamic drag. The system is characterized in the following section.

#### *2.11. Mode Measurement System*

The applied scanning laser vibrometer system (Figure 2) consists of a Polytec OFV-5000 vibrometer controller and Polytec OFV-505 sensor head, two-axis ThorLabs GVS012 scanning mirrors, a National Instruments Compact DAQ-9174 chassis with an NI 9223 analog input card (to ensure the acquisition of the measured and the reference vibrations), and an NI 9263 voltage output card (to adjust the mirrors' positions).

**Figure 2.** The scanning vibrometer setup to used verify that the first resonance mode has been reached.

The system is controlled by a custom-built, PC-based LabVIEW application, which defines and corrects the movement of the mirrors and the laser focus to measure at selected points on the PEH. In addition to these operations, the program acquires the vibrometer velocity data for every set point. The reference vibration was measured with a B&K 4507-B-001 accelerometer to enable the phase compensation of the discrete vibration measurements at the points selected on the PEH's surface. This step ensures the determination of the vibration mode and allows the measurement of the acoustic losses. As the laser vibrometer acquires the velocity of every measured point, it was integrated into the spectral domain to yield the displacement amplitude. The obtained data can be used to confirm the resonance mode and to calculate the acoustic power produced by the PEH.

#### *2.12. Efficiency Measurement System*

The system comprises two parts (Figure 3), one measuring the harvesters' parameters and the other controlling the vibrations. By definition, the measuring component sets the resistance via an Agilent 34970A data acquisition/switch unit and communicates the vibration parameters, frequency, and amplitude to its counterpart. Here, the vibrations are monitored and adjusted by a PI controller, and after they have settled down (error < 0.1% for 10 s, the system will measure the harvester's output voltage on the pre-set load, the displacement of the harvester's tip, and the shaker's current and voltage. The acceleration data are then transferred from the vibration control unit. The smaller displacement amplitudes are measured by a Polytec PDV 100 vibrometer, while in the higher ones a Micro-Epsilon optoNCDT 1401-20 optical sensor is employed. The measurement is performed with a 24-bit, ±5 V, NI 9234 analog input card or, if the measured harvester's output voltage reaches above 5 V, a 24-bit, ±30 V NI 9232 card. To minimize the noise (<40 μV and 127 μV, respectively) in the experiment, we selected sampling frequencies exhibiting the maximum decimation rate of 256. Moreover, the noise and possible scalloping loss were further reduced by extracting the effective value from the signal spectrum via an FFT with a flat top window [24].

**Figure 3.** The system used to measure the power and efficiency.

The vibration control is ensured by a Compact RIO 9067 device equipped with NI 9234 analog input and NI 9263 analog output cards. The generated signal is transmitted to a B&K 2719 power amplifier, which drives a TIRA 52110 vibration shaker. The vibrations are measured by the B&K 4507-B-001 accelerometer, and the signal returns to the Compact RIO via an MFF M28 supply module. The accuracies of the accelerometer and the MMF supply module correspond to 0.4% (in the measured range of 20–300 Hz) and 0.5%, respectively; both of the units are calibrated. To reduce the noise, we employed the method characterized above.

The resistive load is configured by a set of different parallel combinations of resistors on an Agilent 34904A matrix switch inserted into the Agilent 34970A measurement unit. Each of the two switches contains nine resistors, whose parallel combinations cover two decades of resistance (meaning that four decades in total are covered). Throughout the measurement process, the resistances are controlled by an Agilent 34401A multimeter, with the repeatability of the measurements amounting to less than 0.04%. When the parallel combination of the load and the measurement card input resistances produces a significant error, a voltage follower is employed. The entire system is located in an air-conditioned laboratory, where all the measurements take place.

#### *2.13. Comparison with the Methodology Currently in Use*

To facilitate comparison of the proposed efficiency measurement technique with the existing methodology, we calculated the efficiency values in the PEHs at the devices' resonant frequencies and optimum loads by using the standard approach. Even though novel computing formulas are available, such as that described in [25], these share the same simplified principle, which interprets the harvester's mechanical power as the input power. Thus, we decided to compare our new metrics with a well-established and widely used option (Equation (1)) [7] (simplified calculated). The procedure embodied in Equation (1) was experimentally verified by, for instance, the authors of reference [6], whose efficiency measurement concept is also used herein for comparison (simplified measured).

The data enabling the calculation of the aforementioned efficiencies were gathered during merely one measurement, meaning that the error that arises from diverse conditions cannot occur.

#### **3. Results and Discussion**

The measurement system was tested on two commercially available PEHs, one being a PZT PEH MIDE PPA-1011 (PZT PEH), which has a rigid structure with a FR4 base, and the other embodied in a PVDF TE LDTM-028K (PVDF PEH), which exhibits a flexible foil structure with an integrated 720 mg tip mass. The dimensions of these PEHs are shown in Figure 4; the capacitances equal 90.5 nF and 485 pF, respectively.

**Figure 4.** The tested PEHs: the MIDE PPA-1011 (**a**) and the TE LDTM-028K (**b**).

Initially, each resonance frequency was found via the sine sweep, and then the vibration mode (Figure 5) of the examined harvester was checked. The measured amplitudes enabled us to calculate the acoustic power corresponding to the aerodynamic drag. At the final stage, we carried out the same measurement cycles at different vibration amplitudes. The 3D proposed efficiency plots (side and top views) relating to selected vibration amplitudes are shown in Figure 6.

**Figure 5.** The confirmed oscillation: the first longitudinal mode in the MIDE PPA-1011 at 1 grms (**a**) and 2 grms (**b**); an in the TE LDTM-028K at 1 grms (**c**) and 2 grms (**d**).

**Figure 6.** *Cont.*

**Figure 6.** The proposed efficiency rates measured near the harvesters' resonance frequencies, at 1 grms and 2 grms: the zirconate titanate (PZT) PEH, (**a**,**b**, respectively); the polyvinylidene fluoride (PVDF) PEH, (**c**,**d**, respectively). The 3D plots and the top views are shown on the left-hand and the right-hand sides, respectively. Note the different efficiency scales for the two PEHs.

By common definition, the power generated by the PEHs increases with the vibration amplitude. In the PZT PEH, the power varied between 165 μW at 0.25 grms and 3.78 mW at 2 grms. The PVDF PEH, however, is smaller and exhibits a higher optimal load, and thus the power values dropped to between 11.6 nW at 0.25 grms and 793 nW at 2 grms.

#### *3.1. Uncertainties of the Measurement System*

The type B uncertainty was calculated analytically from the devices' specifications. In this uncertainty, the cause consisted in the errors of the NI measurement cards and the accelerometer chain. The error distribution in the individual devices was considered normal (Gaussian), allowing us to calculate the uncertainty. The type B uncertainty varied from 0.003 to 0.4%A for the interval of 1*σ*. The unit of efficiency is usually %, and the same unit normally expresses relative error. In this paper, for clarity, the symbol %A thus denotes the measured efficiency or its absolute error, whereas %R stands for the relative error.

The type A uncertainty was determined as the worst case scenario from all measurements at different frequencies and vibration amplitudes, with each of these uncertainties established from 300 repetitive measuring cycles. The result to express the uncertainty was computed from the post-processed values. In general terms, this method involves the random measurement errors and their filtering by means of an FFT with a flat top window. The relative type A uncertainty amounted to 0.05, 0.03, and 0.002%R in the current, acceleration, and phase, respectively. The absolute combined uncertainty relating to the proposed efficiency measurement ranged between 0.007%A and 0.5%A for the interval of 1*σ*.

#### *3.2. Uncertainty in Simplified Efficiencies*

In the simplified measured and proposed efficiencies, the uncertainties were calculated similarly, from the measurement system's uncertainties. The absolute combined uncertainty relating to the simplified measured efficiency varied between 0.015%A and 0.27%A for the interval of 1*σ*.

The uncertainty of the simplified calculated efficiency, however, arises from the uncertainty of the parameters *k*sys and *Q*, which are calculated from relevant frequencies. For instance, in the PZT PEH, vibrating at 1 grms, the open natural frequency attained *f*<sup>o</sup> = 123.33 Hz and the short one equaled *f*<sup>s</sup> = 123.00 Hz, both at the relative uncertainty 0.12%R. The two frequencies being close to each other, the relative uncertainty of *k*sys amounted to 63.3%R. This rate, above all, then caused the combined absolute uncertainty of the simplified calculated efficiency to be much higher than those in the measured efficiencies, namely, it ranged between 3.05 and 5.56%A for the interval of 1*σ*.

#### *3.3. Efficiency of the PZT Ceramic Based Harvester*

The diagram characterizing the PZT PEH (Figure 7) indicates that the simplified measured and simplified calculated efficiencies yielded similar results; the measured one, however, has a markedly lower uncertainty. These efficiencies are three to five times greater than the proposed efficiency because they did not take into account the mechanical losses.

Moreover, all types of efficiency exhibited a decreasing trend when the vibration amplitude rose. The proposed class nevertheless declined slightly faster due to the growth of the mechanical losses, which occurs with increasing vibration amplitude.

In the given context, Figure 8 shows the individual powers of the PEH. It is worth noting that the aerodynamic loss was higher than the PEH's mechanical power, the difference being approximately 50%. Furthermore, the amount of the vibration input power (*P*VIB) transferred to the mechanical power of the PEH (*P*M) decreased with increasing vibration amplitude from approx. 40% to 20%. This indicates that the mechanical losses increased with rising vibration amplitude.

**Figure 7.** The PZT PEH: Comparing the efficiency values established via our method (proposed) with those calculated from the PEH's parameters for the optimum condition (simplified calculated) and measured from the PEH's mechanical power (simplified measured). The error bars represent the 1*σ* interval.

**Figure 8.** The PZT PEH: Comparing the power delivered by the shaker to the PEH, with the power dissipated in the air, the PEH mechanical vibration power, and the electrical power output.

#### *3.4. Efficiency of the PVDF Foil Based Harvester*

In the PVDF PEH, the simplified calculated and simplified measured efficiencies (Figure 9) amounted to approximately 1.5 %A; the proposed efficiency was significantly smaller, varying between 0.005 and 0.015 %A. Such a condition is caused mainly by the harvester's low stiffness, an effect stemming from the flexible structure. Thus, even though the displacement of the PEH's tip is notable, the vibration power remained very small, being less than 1%R of the input vibration power (Figure 10). Consequently, due the major displacement, the aerodynamic losses produced more than a half of the overall ones. Moreover, the rather high optimum load resistance (330 kΩ) means that the power output was small, ranging between 11.6 nW at 0.25 grms and 793 nW at 2 grms.

**Figure 9.** The PVDF PEH: Comparing the efficiency values established via our method (proposed) (**a**) with those calculated from the PEH's parameters for the optimum condition (simplified calculation) and measured from the PEH mechanical power (simplified measurement) (**b**). The error bars represent the 1*σ* interval.

**Figure 10.** The PVDF PEH: Comparing the power delivered by the shaker to the PEH with the power dissipated in the air, the PEH mechanical vibration power, and the electrical power output.

#### *3.5. Comparing the Efficiencies*

The simplified efficiencies (measured and calculated) exhibited the same results, similarly to the findings presented by Yang et al. [6].

The proposed efficiency showed lower values because the mechanical losses not only need to be taken into account but are also, as is obvious from the measured data, higher than the actual mechanical power of the PEH (*P*M). In general terms, it is then possible to claim that the mechanical losses embody a significant factor affecting the performance of PEHs. Ignoring the losses during in the input power calculations involved in the computation of the simplified efficiency computations of the PEH can cause a significant error when different types of PEHs are compared for various purposes. Our novel method enables direct efficiency measurement and does not omit the mechanical losses of the PEH.

Even in situations where the mechanical losses can be neglected, the calculation of the efficiency from the parameters (1) is affected by a large uncertainty. This arises mainly from the uncertainty of the system coupling coefficient (*k*sys). If the same efficiency is measured via the PEH's mechanical power (*P*M), the uncertainty is significantly lower.

Moreover, excluding the harvester's parameters allows us to use the proposed method for comparing different principles of power generation from ambient vibrations. Our system is capable of measuring a wide range of efficiencies, from tens of percent to units of parts per million, and the shape of the data (Figure 6) appears to be credible even at low efficiency rates. The uncertainty of the measured efficiency is also provided.

#### *3.6. Power Flow*

Once the PEH has been characterized, we can determine the individual power components according to the classification in Figure 1. The individual power components in the absolute value and their contribution to the total power consumed by the PEH are shown in Table 1. These items of information allow us to localize the most significant losses in the power conversion.


**Table 1.** Comparing the individual power components in the tested PEHs at 1 grms.

The most significant types of loss were the acoustic and thermal ones, which together dissipated more than two thirds of the total power input. Regrettably, we were able to find only a few papers that focused on reducing such losses [26–28]. By contrast, many research teams are developing new piezomaterials with the aim to reduce the internal losses; these processes, however, dissipate only a small portion of the power input. Finally, the electrical losses are negligible when compared to the other categories. The intensive research and development of power conditioning circuits will nevertheless allow the electrical losses to be reduced even further in the future.

In conclusion, let us emphasize again that multiple efforts are being made to minimize the internal and electrical losses (which embody only a small portion of the total losses), whereas the type of loss that exhibits the maximum contribution to the total sum remains mostly ignored.

#### **4. Conclusions**

By reviewing the papers that discuss techniques for calculating the efficiency of PEHs, we concluded that researchers regularly compute the input power as the power stored in the harvester's movement. Such a procedure, however, causes the mechanical losses induced by the actual motion to be ignored. We designed a method to measure the efficiency of energy harvesters in such a manner that the mechanical losses generated by internal friction and aerodynamic drag are taken into account. Unlike the traditional approach, our technique relies on computing the input power out of the vibration power extracted from the shaker, with the losses preserved. In the experiments, a mode measurement system was employed to estimate the losses linked to the aerodynamic drag.

We developed a system capable of characterizing the power and measuring directly the overall efficiency of a harvester; the combined absolute uncertainty amounted to less than 0.5 %A. The measurement system was tested on commercially available PZT and PVDF PEHs: a MIDE PPA-1011 and a TE LDTM-028K. In the former, the efficiency at the resonant frequency and optimal load varied from 6.1 %A at 0.25 grms to 2.2 %A at 2 grms. In the latter, the same quantity reached approximately 0.007 %A in all of the measured vibration amplitudes (0.25–2 grms).

The efficiency data relating to the optimum load and frequency were compared with the calculations based on the PEHs' coupling coefficients and quality factors. In the PZT harvester, the proposed measured efficiency rate was roughly three to five times lower than the efficiencies determined via the existing methodology. The PVDF device, however, exhibited rates markedly lower than those established in the PZT one. Generally, it is then possible to claim that the mechanical losses embody a significant factor affecting the performance of PEHs.

We localized the main sources of loss, identifying these with the thermal and acoustic losses caused by the movement of the PEH. The internal losses, arising from the piezoelectric efficiency, consist of only a small portion of the total dissipated power, and the electrical ones, generated by the internal resistance, are negligible. Regrettably, while efforts are being made to minimize those sources of loss that exert a minimum impact on the total losses, only a few research groups focus on the major sources.

As our metrics are independent of the PEHs' parameters, they can be applied in comparing the efficiencies of diverse vibration harvesters. The proposed method may accelerate the development of new vibration harvesters in general because the partial procedures comprised within the technique allow easy and accurate measurement of the devices' performances. By extension, as suggested above, an interesting benefit lies in the possibility of comparing the outcomes delivered by different research teams.

In terms of future research, we intend to improve the measurement system to better measure and model the individual types of loss. Such an aim should then produce an increase in the efficiency of piezo energy harvesters.

**Author Contributions:** Conceptualization, J.K. (Jan Kunz), J.F., S.P. and P.B.; methodology, J.K. (Jan Kunz), S.P. and P.B.; software, J.K. (Jan Kunz) and S.P.; validation, J.F., P.B., J.K. (Jakub Krejci) and Z.H.; formal analysis, J.K. (Jan Kunz) and S.P.; investigation, J.K. (Jan Kunz), J.F., S.P. and P.B.; resources, J.F., S.P., P.B. and Z.H.; data curation, J.K. (Jan Kunz), S.P. and J.K. (Jakub Krejci); writing—original draft preparation, J.K. (Jan Kunz) and S.P.; writing—review and editing, J.K. (Jan Kunz), S.P., J.F., P.B. and Z.H.; visualization, J.K. (Jan Kunz), J.F., S.P. and J.K. (Jakub Krejci); supervision, P.B., S.K. and Z.H.; project administration, S.K. and Z.H.; funding acquisition, S.K. and Z.H. All authors have read and agreed to the published version of the manuscript.

**Funding:** The completion of this paper was made possible by the grant No. FEKT-S-20-6205— "Research in Automation, Cybernetics and Artificial Intelligence within Industry 4.0" financially supported by the Internal science fund of Brno University of Technology. The work was supported by the infrastructure of RICAIP that has received funding from the European Union's Horizon 2020 research and innovation programme under grant agreement No. 857306 and from Ministry of Education, Youth and Sports under OP RDE grant agreement No CZ.02.1.01/0.0/0.0/17\_043/0010085.

**Institutional Review Board Statement:** Not applicable.

**Informed Consent Statement:** Not applicable.

**Data Availability Statement:** The data used for the manuscript are available for researchers on request.

**Acknowledgments:** The authors gratefully acknowledge the contribution of Premysl Dohnal who edited and proofread this paper.

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

#### **Abbreviations**

The following abbreviations are used in this manuscript:


#### **References**


## *Article* **Load Resistance Optimization of Bi-Stable Electromagnetic Energy Harvester Based on Harmonic Balance**

**Sungryong Bae <sup>1</sup> and Pilkee Kim 2,3,\***


**Abstract:** In this study, a semi-analytic approach to optimizing the external load resistance of a bi-stable electromagnetic energy harvester is presented based on the harmonic balance method. The harmonic balance analyses for the primary harmonic (period-1*T*) and two subharmonic (period-3*T* and 5*T*) interwell motions of the energy harvester are performed with the Fourier series solutions of the individual motions determined by spectral analyses. For each motion, an optimization problem for maximizing the output power of the energy harvester is formulated based on the harmonic balance solutions and then solved to estimate the optimal external load resistance. The results of a parametric study show that the optimal load resistance significantly depends on the inductive reactance and internal resistance of a solenoid coil—-the higher the oscillation frequency of an interwell motion (or the larger the inductance of the coil) is, the larger the optimal load resistance. In particular, when the frequency of the ambient vibration source is relatively high, the non-linear dynamic characteristics of an interwell motion should be considered in the optimization process of the electromagnetic energy harvester. Compared with conventional resistance-matching techniques, the proposed semi-analytic approach could provide a more accurate estimation of the external load resistance.

**Keywords:** electromagnetic energy harvester; bi-stable oscillator; load resistance optimization; frequency response analysis; harmonic balance method

#### **1. Introduction**

Energy harvesting is a series of processes that convert ambient energy (vibration, sunlight, wind, raindrop, earth heat, human motion, etc.), usually unused and wasted, into electric energy that is stored in a capacitor or used to supply power to micro-devices. The energy harvesting technology is considered a promising alternative to chemical batteries for autonomous wireless microdevices. In particular, vibration-based energy harvesters have been intensively studied and developed because vibration energy sources exist almost everywhere in our daily life and industrial environment, and its conversion mechanism is relatively simple.

The early design of vibration energy harvesters was based on a linear resonator coupled with electromagnetic, electrostatic, or piezoelectric transducers [1,2]. Such linearresonant energy harvesters produce high electric power only in a narrow resonant frequency band, which is not suitable for harvesting energy from real-world vibration sources with time-varying and/or multiple frequency components [3,4]. Accordingly, energy harvesting systems of various structures (e.g., frequency-tunable oscillators [5–7], an array of multimodal oscillators [8–10], and non-linear mono-stable [11–14] or multi-stable oscillators [15–26]) have been investigated to deal with the narrow bandwidth problem of the linear system.

**Citation:** Bae, S.; Kim, P. Load Resistance Optimization of Bi-Stable Electromagnetic Energy Harvester Based on Harmonic Balance. *Sensors* **2021**, *21*, 1505. https://doi.org/ 10.3390/s21041505

Academic Editors: Zdenek Hadas, Saša Zelenika and Vikram Pakrashi

Received: 10 December 2020 Accepted: 16 February 2021 Published: 22 February 2021

**Publisher's Note:** MDPI stays neutral with regard to jurisdictional claims in published maps and institutional affiliations.

**Copyright:** © 2021 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (https:// creativecommons.org/licenses/by/ 4.0/).

A bi-stable oscillator that possesses a double-well potential energy function is known to be one of the most advantageous candidates to overcome the limitation of the linear system. The dynamic behavior of the bi-stable oscillator is generally more complicated than that of mono-stable oscillators owing to the existence of a saddle barrier between the potential energy wells, but its large-amplitude oscillation (called interwell motion) tends to cover a much broader frequency range. During the past decade, various types of bi-stable energy harvester have been proposed based mainly on magnetically (repulsive/attractive) and/or mechanically buckled structures, and their performance has been evaluated under ambient vibration sources of the swept sine, impulse, and Gaussian white noise types [1,2,18]. The performance of a bi-stable energy harvester depends on the amount of ambient vibration energy transmitted to the energy transducer, that is, transmissibility, which is related to the features of a forced interwell motion, such as the required intensity of excitation, the associated frequency bandwidth, and its oscillation frequency and amplitude. Various optimization strategies have been proposed to enhance the performance of vibration energy harvesters (e.g., artificial intelligence-based optimizations, structural optimizations, and nonlinear optimizations) [27,28]. Additionally, an external load resistance connected to an energy harvester for harvesting electrical energy is also an important design factor that should be optimized for maximizing performance. The resistance-matching technique has been commonly used for estimating the optimal load resistance [29], but it applies in principle to systems where both the source and load devices are linear. For non-linear systems (including a bi-stable system considered in this study), the resistance-matching technique may give an inaccurate estimation of the optimal load resistance with significant error, depending on the system and the excitation conditions.

In this study, a semi-analytic approach based on the harmonic balance method is proposed to optimize the external load resistance of a bi-stable electromagnetic energy harvester accurately, when compared with conventional resistance-matching techniques. To this end, the frequency responses of the primary harmonic and subharmonic interwell motions are numerically obtained first, and their solution forms are assumed in accordance with the results of spectral analyses. For each interwell motion, an optimization problem was formulated and solved, based on the harmonic balance solution of the assumed Fourier series, to calculate the optimal external load resistance of the energy harvester. In addition, a parametric study was conducted to investigate the effects of the system parameters and excitation conditions on the optimal load resistance of the electromagnetic energy harvester, which is not considered in conventional resistance-matching techniques; the results are discussed.

#### **2. Bi-Stable Oscillator Model of an Electromagnetic Energy Harvester**

Figure 1 shows the schematic diagram of an electromagnetic energy harvester considered in this study, which is composed of two parts: a clamped-clamped beam oscillator and an external solenoid coil structure. The inertial mass, along with two permanent magnets, is attached to the middle of the stainless-steel beam to adjust the natural frequency of the beam oscillator. The solenoid coil structure is fixed to the external space and aligned with the inertial mass, so that the permanent magnet can oscillate through the inner hole of the solenoid coil when the base of the energy harvesting system is excited by a certain ambient vibration. The coil is connected to an external load resistance for harvesting electrical energy. With this structure for electromagnetic induction, the energy harvesting system can convert ambient vibration energy to electrical energy that would be used (or stored in a capacitor) for powering micro-devices.

The performance of the energy harvesting system depends significantly on the transmissibility of the system, that is, the amount of ambient vibration transmitted through the beam oscillator to the permanent magnet traveling in the inner hole of the coil. Thus, the electromagnetic energy harvester is designed based on a statically bi-stable oscillator to better respond to various types of ambient vibration with wide-band frequency components. To this end, an axial load is imposed on the beam in the longitudinal direction by

moving the right clamped end slightly to the left, which causes the buckling phenomenon of the beam. It is well known that at the moment when buckling occurs, new symmetric non-trivial equilibria (*x* = ±*x*0) start bifurcating from the trivial one (*x* = 0), which is called the pitch-fork bifurcation of equilibrium. After the buckling, the beam structure becomes unstable at the trivial equilibrium, and is stable at the two non-trivial upper and lower equilibria depicted in Figure 1.

**Figure 1.** Schematic diagram of an electromagnetic bi-stable energy harvester composed of a clampedclamped beam and an external solenoid coil. The right end of the beam is movable in the horizontal direction but can also be fastened at a certain position. The gravitational field is assumed to be in the perpendicular direction to the given plane.

The mathematical model for the aforementioned energy harvesting system is given by the following electro-magneto-mechanical oscillator model [30–32]:

$$
\ddot{\mathbf{x}} + 2\zeta\omega\_0 \dot{\mathbf{x}} - \gamma \mathbf{x}\_0^2 \mathbf{x} + \gamma \mathbf{x}^3 + \kappa\_0 I = -A\_b \cos \Omega t,\tag{1a}
$$

$$L\_0 \dot{I} + (R\_{\text{f}} + r\_{\text{L}})I = \beta \dot{\mathbf{x}}\_{\text{\prime}} \tag{1b}$$

where *x* is the displacement of the inertial mass attached to the middle of the beam; *x*<sup>0</sup> indicates the distance of the non-trivial equilibria measured from the trivial one; *γ* is the coefficient given by (*ω*0/4*x*0) 2 ; *I* is the electrical current flowing across the load resistance *Re*; *ω*<sup>0</sup> and *ζ* are the natural angular frequency and damping ratio, respectively; *κ*<sup>0</sup> and *β* are the electro-magneto-mechanical coupling constants; *rL* and *L*<sup>0</sup> are the internal resistance and inductance of the coil, respectively; and *Ab* and Ω are the magnitude and frequency of the base acceleration, respectively. The baseline parameter values of the oscillator model are listed in Table 1. In addition, for the purpose of generality, the following dimensionless variables and parameters are introduced to Equations (1a) and (1b):

$$\text{If } y = \frac{\text{x}}{\text{x}\_0 \text{\textdegree I}}, \text{ I} = \frac{\text{I}}{\text{I}\_0}, \text{ } \pi = \omega\_0 t, \text{ } A = \frac{A\_b}{\text{x}\_0 \omega\_0^2 \text{\textdegree I}}, \omega\_\Gamma = \frac{\omega}{\omega\_0}, \text{ } \pi = \frac{\kappa\_0 \beta L\_0}{a^2 r\_L^2}, \text{ } a = \frac{\omega\_0 L\_0}{r\_L}, \text{ } R = \frac{R\_\varepsilon}{r\_L}, \tag{2}$$

and the resulting dimensionless forms of governing equations can be obtained in the form:

$$
\ddot{y} + 2\zeta \dot{y} - \frac{1}{2}y + \frac{1}{2}y^3 + \kappa f = -A\cos\omega \tau,\tag{3a}
$$

$$
\alpha \dot{f} + (\mathcal{R} + 1)f = \dot{g}\_{\prime} \tag{3b}
$$

where *ω* are and the dimensionless forcing frequency (or the frequency ratio); *α* is the dimensionless resonant inductive reactance (or the ratio of the resonant inductive reactance of the coil, *ω*0*L*0, to the internal resistance of the coil, *rL*); and *R* is the dimensionless external load resistance (or the ratio of the external load resistance to the internal resistance of the coil).


**Table 1.** Parameter values of the electro-magneto-mechanical oscillator [30,32].

#### **3. Numerical Frequency Response Analysis of an Electro-Magneto-Mechanical Oscillator Model**

#### *3.1. Steady-State Response and Its Bifurcation*

In this work, a series of numerical simulations of the electro-magneto-mechanical oscillator model expressed in Equations (3a) and (3b) are performed to investigate the dynamic behaviors of the system and the associated energy harvesting performance. By introducing the state vector [*y*1, *<sup>y</sup>*2, *<sup>y</sup>*3, *<sup>y</sup>*4] <sup>=</sup> . *y*, *y*, *J*, *θ* ! , where *θ* = *ωτ*, into Equations (3a) and (3b), the dimensionless governing differential equations can be rewritten into an autonomous system of first-order differential equations in the following vector form:

$$\dot{\mathbf{y}} = f(\mathbf{y}, t) = \begin{bmatrix} -2\tilde{\mathbf{y}}y\_1 + \frac{1}{2}y\_2 - \frac{1}{2}y\_2^3 - \kappa y\_3 - A\cos y\_4 \\ y\_1 \\ -\frac{R+1}{a}y\_3 + \frac{1}{a}y\_1 \\ \omega \end{bmatrix}. \tag{4}$$

Frequency response analyses are carried out by the direct numerical integration of Equation (4) based on the Runge–Kutta methods. For all the numerical simulations, the baseline parameter values given in Table 1 are used, unless otherwise stated. The steadystate forced responses of the electromagnetic energy harvester are evaluated first, while the frequency of the base excitation varies by means of frequency marching or frequency sweeping, and then their Floquet multipliers are obtained by solving the eigenvalue problem of the following monodromy matrix:

$$\mathbf{M} = \prod\_{i=N}^{1} \exp\left[\int\_{\tau\_{i-1}}^{\tau\_i} \mathbf{J}(\tau) d\tau\right],\tag{5}$$

where **J** is the Jacobian matrix, and *N* is the number of equally spaced time intervals included within one period. The bifurcation type of the stationary oscillation can be determined from the investigation of its Floquet multipliers. When the stability of the steady-state solution changes, the accompanying bifurcations can be classified into three different types: saddle-node bifurcation (if the largest Floquet multiplier *μ* is positive 1, i.e., *μ* = +1), period-doubling bifurcation (*μ* = −1), and Neimark–Sacker bifurcation (|*μ*| = 1 and Im[*μ*] = 0).

#### *3.2. Non-Linear Oscillations Owing to Double-Well Potential: Intrawell and Interwell Motions*

The electromagnetic bi-stable energy harvester possesses a double-well potential energy function, as shown in Figure 2a. There are one trivial saddle point and two non-trivial center points in the potential function. The two potential wells are in symmetric configuration with respect to the potential saddle barrier formed at the trivial point of *y* = 0. This potential saddle barrier leads to two different types of forced response: intrawell and interwell motions, depending on the intensity of the base excitation. When the base excitation is weak, the system only oscillates with small amplitude inside one of the two potential wells (as demonstrated by the thin line in Figure 2b), which is called intrawell motion, owing to the existence of the potential barrier. When the base excitation becomes relatively strong, the system gains enough kinetic energy to cross the potential barrier (the thick line

in Figure 2b). Accordingly, it begins oscillating with a much larger amplitude, producing higher electrical power, and it is called interwell motion. The interwell motion tends to occur over a much broader frequency band than the resonant frequency band of linear oscillators. Using this interwell motion, the electromagnetic bi-stable energy harvester can harvest energy from ambient vibration energy sources, including wide-band frequency components. The aforementioned high-energy interwell motion of the electromagnetic energy harvester appears in various types of oscillations: primary harmonic motion, its subharmonic and superharmonic motions, and even non-periodic chaotic motion. Among these interwell motions, the primary harmonic motion of period *T* = 2*π*/*ω* has been used the most for broadband energy harvesting applications. However, it was recently reported that subharmonic motions are also useful and suitable [29].

**Figure 2.** (**a**) Double-well potential energy function of the electromagnetic bi-stable energy harvester. The open and solid circles indicate the saddle and center points, respectively. (**b**) Phase plane representation of the intrawell and interwell oscillations (of period *T* = 2*π*/*ω*), which are denoted by thick (blue) and thin (magenta) lines, respectively, obtained when *A* = 0.02 and *ω* = 0.82. In (**b**), the red points indicate the stroboscopic projections of the oscillations, which are synchronized with the period *T* of the base excitation.

Figure 3 shows the frequency responses for (a) displacement and (b) output power of the electromagnetic energy harvester obtained when the base acceleration is 0.02. In this figure, the subharmonic interwell motions of periods 3*T* and 5*T* and the primary harmonic interwell motion of period *T* are presented along with the intrawell motion for comparison purposes. The subharmonic interwell motions are much larger in amplitude than the intrawell motion (as shown in Figure 3a), although they are smaller than the primary harmonic interwell motion, and accordingly can produce relatively higher output powers (as shown in Figure 3b). This is more clearly demonstrated in Figure 4 by directly comparing the period-1*T*, 3*T*, and 5*T* interwell motions on the state plane. The period-3*T* interwell motion has a notably wide frequency band, similar to that of the period-*T* interwell motion. It is interesting that the frequency bands of the period-1*T*, 3*T*, and 5*T* interwell motions are continuously connected to each other and tend to form a whole wide frequency band. The above observations indicate that the subharmonic interwell motions in addition to the primary harmonic interwell motion can play an important role in enhancing the operating frequency band of electromagnetic energy harvesters.

**Figure 3.** Frequency responses for (**a**) displacement and (**b**) output power of the electromagnetic bi-stable energy harvester obtained when *A* = 0.02: (1*T*) primary harmonic period-*T* interwell motion, (1*T*-SO) period-*T* intrawell motion, (3*T*) subharmonic period-3*T* interwell motion, and (5*T*) period-5*T* interwell motion. The saddle-node bifurcations and perioddoubling bifurcation are designated by points SB and PD, respectively. In (**a**), the red dots indicate the stroboscopic projections of the oscillations, which are synchronized with the period *T* of the base excitation. The dashed-dotted lines are the harmonic balance solutions for the primary harmonic and subharmonic interwell motions, which will be discussed in Section 4.

**Figure 4.** Phase plane representations of (**a**) period-3*T* and (**b**) period-5*T* motions obtained with the base acceleration of *A* = 0.02 and forcing frequencies of (**a**) *ω* = 1.37 and (**b**) *ω* = 1.65, respectively. The period-*T* interwell and intrawell motions of Figure 2b are also presented for comparison purpose. The solid circle dots indicate the stroboscopic projections of the oscillations, which are synchronized with the period *T* of the base excitation.

#### **4. Semi-Analytic Approach to Optimizing the External Resistance**

The external load resistance of an electromagnetic energy harvester, together with the transmissibility, is an important design factor, which must be optimized to maximize its output electrical power. An impedance-matching technique has been commonly used to solve this optimization problem, with the assumption that the inductance of the coil is negligible. In this case, the external load resistance is actually matched with the internal resistance of the solenoid coil; thus, it will be called 'resistance matching' in this study. However, it should be noted that the impedance- (or resistance-) matching technique only applies when both the source and load devices are linear, although they are frequently used for non-linear systems at the cost of accuracy. In this work, a semi-analytical approach based on the method of harmonic balance is implemented to evaluate the optimal value of

the external load resistance accurately, with which the interwell motion of the system can most efficiently produce electrical power.

First, Fourier analyses are performed to find the appropriate solution forms of the steady-state primary harmonic and subharmonic interwell motions, which depend on the type and intensity of non-linearity [33], and the results of the time responses and the corresponding amplitude spectra are presented in Figure 5. The forcing period of the base excitation is given by *T* = 2*π*/*ω*, and accordingly the oscillation frequencies of the subharmonic period-*3T* and 5*T* motions tend to be lower than the forcing frequency *ω*, as shown in Figure 5c,e. Such subharmonic behavior is known to be a typical feature of a bi-stable system. Therefore, the fundamental oscillation frequencies of the interwell motions can be defined by *ω<sup>k</sup>* = *ω*/*k*, where *k* = 1, 3, and 5 for the period-1*T*, *3T*, and 5*T* motions, respectively. From Figure 5b,d,f, it can be observed that the interwell motions only have two or three important frequency components: *ω*1- and 3*ω*1-components for the period-*T* motion, *ω*3- and 3*ω*3-components for the period-*3T* motion, and *ω*5-, 3*ω*5- and 5*ω*5-components for the period-5*T* motion, as shown in Figure 5b,d,e, respectively. Thus, the steady-state periodic solutions can be assumed in the following forms of the *N*th-order Fourier series:

$$\text{(Period } -T \text{ motion) } y = \sum\_{n=1,3} a\_{\text{ll}} \cos n\omega \tau + b\_{\text{ll}} \sin n\omega \tau, \text{ } f = \sum\_{n=1,3} c\_{\text{ll}} \cos n\omega \tau + d\_{\text{ll}} \sin n\omega \tau \tag{6a}$$

$$\text{(Period } -3T \text{ motion) } y = \sum\_{n=1,3} a\_n \cos \frac{n\omega}{3} \pi + b\_n \sin \frac{n\omega}{3} \pi \text{ , } f = \sum\_{n=1,3} c\_n \cos \frac{n\omega}{3} + d\_n \sin \frac{n\omega}{3} \tag{6b}$$

$$\left(\text{Period} - 5T \text{ motion}\right)y = \sum\_{n=1,3,5} a\_n \cos\frac{n\omega}{5} \pi + b\_n \sin\frac{n\omega}{5} \pi \ .\ / = \sum\_{n=1,3,5} c\_n \cos\frac{n\omega}{5} + d\_n \sin\frac{n\omega}{5} \tag{6c}$$

where the highest order of harmonics is given as *N* = 3 for the period-1*T* and 3*T* motions and *N* = 5 for the period-5*T* motion, and all the constant terms are set to zero as the main interest is on the interwell motion, of which the center of oscillation is the trivial point (*y* = 0).

Substituting Equations (6a)–(6c) into Equations (3a) and (3b), followed by equating the coefficients of like harmonic components, leads to a system of non-linear algebraic equations for unknown variables *an*, *bn*, *cn*, and *dn*, which are numerically solved by the Newton–Raphson method. As shown in Figure 3, the harmonic balance solutions (dasheddotted lines) for primary harmonic and subharmonic interwell motions are well matched with the numerical solutions. Additionally, time-domain comparisons are made in Figure 6, where numerical and analytical results are also in a good agreement with each other.

**Figure 5.** *Cont*.

**Figure 5.** (**first column**) Time responses and (**second column**) the associated amplitude spectra, obtained by the fast Fourier transform. In (**a**,**b**) the period-1*T* motion when *A* = 0.02 and *ω* = 0.82, in (**c**,**d**) the period-3*T* motion when *A* = 0.02 and *ω* = 1.37, and in (**e**,**f**) the period-5*T* motion when *A* = 0.02 and *ω* = 1.65. The base excitation used for this simulation is also presented for comparison. In (**a**,**c**,**e**), the red dots indicate the stroboscopic projections of the oscillations, which are synchronized with the period *T* of the base excitation.

**Figure 6.** *Cont*.

**Figure 6.** Comparisons between numerical and analytical solutions for (**a**) period-1*T*, (**b**) 3*T*, and (**c**) 5*T* interwell motions, which are obtained when the amplitude and frequency of the base acceleration are given as *A* = 0.02 and *ω* = (**a**) 0.82, (**b**) 1.37, and (**c**) 1.65.

Using the obtained harmonic balance solutions, the average output power *Pav* of the electromagnetic energy harvester can be evaluated as *P*<sup>2</sup> *av* = ∑ (*c*<sup>2</sup> *<sup>n</sup>* + *d*<sup>2</sup> *<sup>n</sup>*)/2, and the optimal value of the external load resistance is estimated by solving the optimization problem of *∂Pav*/*∂R* = 0. The optimal external load resistance reduces to 1 for the same problem when the resistance matching technique is used.

Figure 7 shows the average output powers, harvested from (a) period-1*T*, (b) period-3*T*, and (c) period-5*T* motions, with respect to the external load resistance. The numerical and semi-analytical results are compared to each other in Figure 7. As shown in this figure, the analytic and numerical results for the average output powers are well matched with each other. Furthermore, the optimal values of the external load resistance were evaluated, and it was observed that the analytic results are in an excellent agreement with numerical results with a maximum relative error of less than (a) 0.3%, (b) 2.0%, and (c) 1.5%. This observation supports the validity of the semi-analytic approach to optimize the external load resistance implemented in this study.

**Figure 7.** *Cont*.

**Figure 7.** Plots of the external load resistance versus the average output power harvested from (**a**) period-1*T*, (**b**) period-3*T*, and (**c**) period-5*T* motions of the electromagnetic energy harvester. The results are obtained by numerical integration (blue solid lines) and the semi-analytic approach (red dashed lines) and compared to each other. For each motion, the comparisons in between are made with three different forcing conditions: (**a**) (*A*, *ω*) = (0.012, 0.6), (0.02, 0.7), (0.03, 0.8), (**b**) (*A*, *ω*) = (0.012, 1.0), (0.02, 1.2), (0.03, 1.4), (**c**) (*A*, *ω*) = (0.012, 1.4), (0.02, 1.6), (0.03, 1.8). The solid circle points indicate the optimal conditions for the external load resistance with which the output power is maximized.

In general, the subharmonic period-3*T* and 5*T* interwell motions tend to produce relatively lower electric power than the primary harmonic period-1*T* interwell motion (but relatively higher than intrawell motions). However, the frequency bands of the period-1*T*, 3*T*, and 5*T* interwell motions are possibly different from each other (as shown in Figure 3), and each of the interwell motions could be used to harvest energy in the associated frequency band, depending on excitation conditions.

#### **5. Effects of Inductive Reactance and Excitation Conditions on Optimal External Resistance**

Using the semi-analytic method mentioned in the previous section, a parametric study is conducted to investigate the effects of the dimensionless resonant inductive reactance *α* (defined by the ratio of *ω*0*L*<sup>0</sup> to *rL* in Equation (3b)) and the excitation conditions (the frequency and intensity of the base excitation) on the optimal load resistance of the electromagnetic energy harvester, which is not considered in the resistance-matching technique.

Figure 8 shows the optimal values of the external load resistance for (a) period-1*T*, (b) period-3*T*, and (c) period-5*T* motions of the electromagnetic energy harvester, obtained by using the semi-analytic method, with respect to the resonant inductive reactance. In this figure, the results clearly indicate that the optimal load resistance of the electromagnetic energy harvester can be significantly affected by the resonant inductive reactance and excitation frequency. For all the period-1*T*, 3*T*, and 5*T* motions, the optimal load resistance increases as the dimensionless resonant inductive reactance increases—-in other words, as the inductance (*L*0) of the solenoid coil becomes relatively larger—-and accordingly, it deviates from the optimal value (i.e., 1) estimated by the resistance-matching technique, in which the inductance of the coil is assumed to be zero. Such a deviation tends to be larger when the frequency (*ω*) of the base excitation is higher, which is more remarkable for the interwell motion of the shorter oscillation period (in the order of periods (a) 1*T*, (b) 3*T*, and (c) 5*T* in Figure 8) or higher oscillation frequency (in the order of (a) *ω*1, (b) *ω*3, and (c) *ω*5). From these observations, it is deduced that the optimal load resistance depends on the inductive reactance of the coil, which is defined by *ωkL*0, in addition to the internal resistance (*rL*) of the coil.

**Figure 8.** Plots of the resonant inductive reactance versus the optimal load resistance for (**a**) period-1*T*, (**b**) period-3*T*, and (**c**) period-5*T* motions of the electromagnetic energy harvester. For each motion (colored thick lines), the results are calculated by using the semi-analytic method, with three different values of excitation frequency and a constant base acceleration of 0.02, and compared with that obtained by the resistance-matching technique (grey thin line).

Figure 9 shows the optimal values of the external load resistance for the period-1*T*, 3*T*, and 5*T* motions of the electromagnetic energy harvester, obtained with three different values of base acceleration: (a) *A* = 0.012, (b) *A* = 0.02, and (c) *A* = 0.03. The calculations are made only in the range of the excitation frequency in which the interwell (period-1*T*, 3*T*, and 5*T*) motions exist. As shown in this figure, the period-1*T* and period-3*T* motions occur in a broad range of excitation frequencies, whereas the period-5*T* motion occurs in a relatively narrow range, which means that the former are more important design factors in broadband energy harvesting applications. The frequency range for each interwell motion tends to be broader as the intensity of the base excitation increases, slightly shifting to a higher frequency. For the period-1*T* motion, the optimal load resistance obviously increases with excitation frequency, whereas for other interwell motions, only small variations in the optimal load resistance are observed in the overall range of excitation frequency (particularly for period-5*T* motion, the resistance is almost constant). This is because the lower the oscillation frequencies of the subharmonic period-3*T* and 5*T* motions (*ω*<sup>3</sup> = *ω*/3 and *ω*<sup>5</sup> = *ω*/5, respectively) are, the weaker the effects of the inductive reactance of the solenoid coil on the optimal resistance. Furthermore, the optimal load resistance is not sensitive to changes in the intensity of the base excitation. As shown in Figure 10, the optimal load resistance for the period-1*T* and 3*T* motions tends to increase slightly with the base acceleration and, contrarily, the one for the period-5*T* motion decreases slightly. However, these changes in the optimal resistance are not significant (less than 0.01%).

**Figure 9.** Plots of the excitation frequency versus the optimal load resistance for the interwell motions (with periods 1*T*, 3*T*, and 5*T*) of the electromagnetic energy harvester. In this simulation, the acceleration of base excitation is set to (**a**) 0.012, (**b**) 0.02, and (**c**) 0.03.

**Figure 10.** *Cont*.

**Figure 10.** Base acceleration versus optimal load resistance for (**a**) period-1*T*, (**b**) period-3*T*, and (**c**) period-5*T* motions of the bi-stable electromagnetic energy harvester.

#### **6. Broadband Energy Harvesting Applications**

The output power of the bi-stable electromagnetic energy harvester is evaluated, while the frequency of the harmonic base excitation varies by means of frequency marching, and compared with that of its linear counterpart in Figure 11 to demonstrate its broadband energy harvesting performance. As shown in this figure, for the linear system, a sharp resonant peak of the output power response appears in a very narrow frequency band. In this case, the natural frequency of the system must be designed to be synchronized with the forcing frequency of the harmonic excitation, otherwise it is likely to fail to produce high output power. On the contrary, the bi-stable system possesses the multiple branches of the interwell motions (including the period-1*T*, 3*T*, and 5*T* motions) in a broad frequency band and thereby can keep normally operating with high output power, regardless of the synchronization of the frequencies. Therefore, the bi-stable energy harvester is definitely more advantageous in real applications, in which the ambient harmonic vibration possibly varies, than the linear system. Actually, there are various types of ambient vibration source such as harmonic, periodic, impulsive, and random excitations. The external load resistance of the bi-stable electromagnetic energy harvester should be optimized in an appropriate manner that depends on the type of the ambient vibration source.

**Figure 11.** Frequency marching responses for the output powers of the electromagnetic bi-stable energy harvester and its linear counterpart, which are obtained when (**a**) *A* = 0.02 and (**b**) 0.03.

Basically, the semi-analytical optimization approach proposed in this study is suitable for harmonic vibration sources, since it is based on the harmonic balance solution. For all the above simulations, the optimal load resistance was evaluated and investigated, when the frequency of the harmonic base excitation was constant, to demonstrate the dynamic characteristics of the optimal resistance. In this section, the proposed optimization approach is further extended for piecewise harmonic excitations with constant amplitude but stepwise changes in frequency.

The general form of the average output power is given in the form:

$$P\_{av} = f\_{RMS}^2 R = \frac{1}{\Delta T} \int\_{T\_0}^{T\_N} f^2 R \, d\tau,\tag{7}$$

where *JRMS* is the root mean square (RMS) of the output current flowing across the external load resistance over the whole-time interval Δ*T*(= *Tn* − *T*0). For most base excitations, the numerical integration method can be commonly used to evaluate the average power of Equation (7). Particularly for the piecewise harmonic excitation, the frequency of which is defined on a sequence of sub-time intervals <sup>Δ</sup>*Tk* = *Tk* − *Tk*−<sup>1</sup> (*<sup>k</sup>* = 1, 2, ··· , *<sup>N</sup>*), the resulting average output power is rewritten as the weighted sum of *N* RMS values:

$$P\_{d\overline{\nu}} = \frac{\Delta T\_1}{\Delta T} \left( \frac{1}{\Delta T\_1} \int\_{T\_0}^{T\_1} f^2 R d\tau \right) + \frac{\Delta T\_2}{\Delta T} \left( \frac{1}{\Delta T\_2} \int\_{T\_1}^{T\_2} f^2 R d\tau \right) + \dots + \frac{\Delta T\_N}{\Delta T} \left( \frac{1}{\Delta T\_N} \int\_{T\_{N-1}}^{T\_N} f^2 R d\tau \right), \tag{8a}$$

or

$$P\_{\theta\overline{\nu}} = \sum\_{k=1}^{n} r\_k f\_k^2 R \text{ with } r\_k = \frac{\Delta T\_k}{\Delta T} \text{ and } f\_k = \sqrt{\frac{1}{\Delta T\_k} \int\_{T\_{k-1}}^{T\_k} f^2 d\tau\_\gamma} \tag{8b}$$

where *Jk* is the RMS output current on the subinterval Δ*Tk*, and *rk* is the weighting factor. Assuming that the base excitation is harmonic on each subinterval Δ*Tk*, the RMS output current (*Jk*) and the associated average power (*J*<sup>2</sup> *<sup>k</sup> R*) can be evaluated by using the harmonic balance solutions of Equations (6a)–(6c). The optimization problem for maximizing the total output power given by Equation (8) is readily formulated and solved to estimate a single optimal value of the external load resistance.

Figure 12a shows the amplitude of the output current response to (red line) piecewise harmonic excitation, which is obtained numerically for the period-1*T* interwell motion of the bi-stable electromagnetic energy harvester. In this simulation, the forcing frequency of the piecewise base excitation decreased from 1.0 by 0.05 to 0.5 for every 1000 cycles, but the amplitude of the base acceleration is set to be 0.02 g. The resulting output current response tends to intermittently vary with stepwise change in excitation frequency. Additionally, the average output power is evaluated with the variation of the external load resistance, as illustrated in Figure 12b. The optimal load resistance, 1.111, is estimated and compared with the analytical one, 1.104, evaluated by using Equations (6) and (8). The relative error of the analytical optimal resistance is very small (less than 0.7%), which supports the assertion that the proposed optimization approach is also valid for the piecewise harmonic excitation. In Figure 12a,b, the results for (blue line) swept-sine excitation are presented for comparison purposes. The swept-sine excitation continuously decreased with the slow sweep rate of 5 × <sup>10</sup>−<sup>6</sup> in the same frequency range. The output current response to the swept-sine excitation gradually decreases in the similar tendency to the piecewise harmonic case but with the difference in amplitude, which results in the difference in average output power. However, the optimal resistance, 1.139, for the swept-sine excitation (particularly with a slow sweep rate) is very similar to that for the piecewise harmonic case, with a small difference (less than 2.6%). This means that the semi-analytical optimal resistance is possibly used as an approximate for slowly swept excitation.

As another example, Figure 13 shows the output current responses to piecewise harmonic and swept-sine excitations, which is obtained numerically for the period-3*T* interwell motion. The forcing frequency of the piecewise base excitation increased from 1.0 by 0.05 to 1.5 for every 1000 cycles. The frequency of the swept-sine excitation increased with the sweep rate of 8.9 × <sup>10</sup>−6. The numerical optimal resistances for the piecewise harmonic and swept-sine excitations are nearly close to each other, approximately 1.056.

The semi-analytical optimal resistance, 1.059, is in good agreement with the numerical results for both the piecewise harmonic and swept-sine excitations.

**Figure 12.** (**a**) Amplitudes of the output current responses to (red line) piecewise harmonic and (blue line) swept-sine excitations, which are obtained for the period-1*T* interwell motions of the bi-stable electromagnetic energy harvester, and (**b**) average output powers evaluated with the variation of the external load resistance.

**Figure 13.** (**a**) Amplitudes of the output current responses to (red line) piecewise harmonic and (blue line) swept-sine excitations, which are obtained for the period-3*T* interwell motions of the bi-stable electromagnetic energy harvester, and (**b**) average output powers evaluated with the variation of the external load resistance.

#### **7. Conclusions**

In this study, a semi-analytic approach to optimizing the external load resistance of an electromagnetic energy harvester was proposed for primary harmonic (period-1*T*) and subharmonic (period-3*T* and 5*T*) interwell motions. The harmonic balance solutions of three interwell motions were obtained first and used to formulate an optimization problem for the external load resistance. The optimal load resistance was obtained, by solving the formulated problem, and investigated. A parametric study was performed to investigate the effects of the system parameters and excitation conditions on the optimal load resistance of the electromagnetic energy harvester. The results showed that for all the period-1*T*, 3*T*, and 5*T* motions, the optimal load resistance tended to increase as the inductance of a solenoid coil was larger and as the frequency of the base excitation was higher, but it was not sensitive to the intensity of the base excitation. Additionally, the optimal load resistance significantly depended on the oscillation frequency (or period) of

the interwell motion—the higher the oscillation frequency (or the shorter the period) is, the larger the optimal resistance; thus, for the period-1*T*, 3*T*, and 5*T* motions in order. All of the above observations consistently indicate that the effect of the inductive reactance of the solenoid coil on the optimal load resistance becomes significant, particularly when the frequency of ambient vibration is relatively high. In this case, the non-linear dynamic behaviors of the interwell motions and the associated inductive reactance of the solenoid coil should be considered as important design factors in the optimization process of the electromagnetic energy harvester. Therefore, when compared with conventional resistancematching techniques, the semi-analytic approach presented in this study could provide a more accurate estimation of the optimal load resistance.

**Author Contributions:** Conceptualization, S.B. and P.K.; methodology, P.K.; software, S.B. and P.K.; validation, S.B.; formal analysis, P.K.; writing—original draft preparation, S.B. and P.K.; supervision, P.K.; project administration, P.K. All authors have read and agreed to the published version of the manuscript.

**Funding:** This work was supported by the National Research Foundation of Korea (NRF) grant funded by the Korean Government (MSIT) (No. NRF-2019R1C1C1009732). This research was also supported by "Research Base Construction Fund Support Program" funded by Jeonbuk National University in 2017.

**Institutional Review Board Statement:** Not applicable.

**Informed Consent Statement:** Not applicable.

**Data Availability Statement:** Data sharing not applicable.

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

#### **References**


#### *Article*
