*Article* **Dynamic Modeling and Experimental Validation of an Impact-Driven Piezoelectric Energy Harvester in Magnetic Field**

#### **Chung-De Chen \*, Yu-Hsuan Wu and Po-Wen Su**

Department of Mechanical Engineering, National Cheng Kung University, Tainan City 701, Taiwan; n16064593@gs.ncku.edu.tw (Y.-H.W.); n16084462@gs.ncku.edu.tw (P.-W.S.)

**\*** Correspondence: cdchen@mail.ncku.edu.tw

Received: 29 September 2020; Accepted: 26 October 2020; Published: 29 October 2020

**Abstract:** In this study, an impact-driven piezoelectric energy harvester (PEH) in magnetic field is presented. The PEH consists of a piezoelectric cantilever beam and plural magnets. At its initial status, the beam tip magnet is attracted by a second magnet. The second magnet is moved away by hand and then the beam tip magnet moves to a third magnet by the guidance of the magnetic fields. The impact occurs when the beam motion is stopped by the third magnet. The impact between magnets produces an impact energy and causes a transient beam vibration. The electric energy is generated by the piezoelectric effect. Based on the energy principle, a multi-DOF (multi-degree of freedom) mathematical model was developed to calculate the displacements, velocities, and voltage outputs of the PEH. A prototype of the PEH was fabricated. The voltages outputs of the beam were monitored by an oscilloscope. The maximum generated energy was about 0.4045 mJ for a single impact. A comparison between numerical and experimental results was presented in detail. It showed that the predictions based on the model agree with the experimental measurements. The PEH was connected to a diode bridge rectifier and a storage capacitor. The charges generated by the piezoelectric beam were stored in the capacitor by ten impacts. The experiments showed that the energy stored in the capacitor can light up the LED.

**Keywords:** piezoelectric energy harvesting; vibration; frequency-up conversion

#### **1. Introduction**

The piezoelectric energy harvesters (PEHs) have received attention in the past two decays due to the high cost for battery replacements in the wireless sensor networks [1]. In the early development of PEHs, the resonators, such as cantilever beams, are used to harvest vibration energy from the environments [2–4]. The narrow bandwidth limits the application of the resonance-based energy harvesters. Many researchers have developed broad-band strategies to broaden the bandwidth of the resonance-based energy harvester, such as nonlinear oscillators, arrayed oscillators and multi-mode coupled oscillators [5].

In recent years, harvesting energy from human motions has become a hot topic because of the increasing demands of wearable devices [6]. The frequency of the human motion is extreme low, typically less than 5 Hz. For such low frequency, the energy conversion efficiency is also low for resonance-based energy harvesters. Although some broad-band strategies are adopted to enhance the efficiency, it is difficult to widen the bandwidth covering a frequency less than 5 Hz [7]. As an alternative, the concept of frequency-up conversion has been proven to be able to harvest energy at high efficiency from energy sources with extreme low frequencies. Umeda et al. presented a concept for energy harvesting based on the frequency up-conversion [8]. They investigate a ball free falling onto a piezoelectric beam to produce a transient vibration. Gu and Livemore presented a two-beam

assembly [9]. The long beam vibrates in responding to the environment with a low frequency and hits the short piezoelectric beam with a high natural frequency. The electrical energy can be extracted from the transient vibration of the short beam with high efficiency. In recent years, some similar studies in energy harvesting based on frequency-up conversion with impact type can be also seen in [10–13].

In the above mentioned frequency up-conversion energy harvesters, the conversions from low frequency to high frequency through mechanical contact. Wickenheiser and Garcia presented a cantilever beam with magnet plucked by magnetic forces [14]. Although it avoids some problems such as surface wear and noises due to contact, the energy generated by the magnetic plucking is usually smaller than those based on mechanical contact. The frequency up-conversion driven by magnetic plucking forces was also adopted in some papers [15,16].

The frequency up-conversion concept has been proven having potential to harvest energy with extreme low frequency, and hence can be adopted to harvest energy from human motion. Pozzi and Zhu presented a device mounted at the knee to generate electrical energy from knee motion [17]. Their device is composed of a stator with plectra and a rotor with cantilever piezoelectric beams. As the knee moves, the stator and rotor have a relative rotational motion and the plectra pluck the beams. A refined device to harvest energy from knee motion was proposed by Kuang et al. [18]. The device is similar to that in [17] except that the beams were plucked by magnetic force. Their experiments showed that a maximum power of 4.5 mW can be harvested from walking. Wei et al. proposed an impact type energy harvester [19]. The presented this device mounted on a human leg to harvest energy from human walk. Their experiments showed that a power of 51 μW was generated for a walking speed of 5 km/h.

The magnetic forces are commonly used in the PEH designs. In the early development, the magnetic forces were assumed as simple forms such as inverse-square model [20–22]. These simple models are easy to implement. However, these forms valid only for some specific displacement range. For the case that the distance between two magnets are short or long, these simple models give inaccurate predictions. Some researchers assumed that the magnet as dipole without occupying a volume [23–26]. For the case of the large distance between two magnets, the model gives an accurate magnetic force. However, when one magnet approaches another, the assumption of the point-dipole become doubtful. Some 3D models have been proposed to calculate the interaction force between two magnets [27,28]. Comparing to the assumed form model and the point-dipole model, the 3D model requires more computation resources as it gives a more accurate prediction in magnetic force. In recent years, the 3D model has been introduced in developing the mathematical models of PEH [16,29–31].

In this paper, an impact-driven piezoelectric energy harvester (PEH) in the magnetic field is presented. A multi-DOF mathematical is developed to investigate the dynamic behaviors of the PEH. A 3-D magnetic force model is also introduced to calculate the magnetic forces between magnets. The voltage responses and energy harvested by the PEH can be calculated by the model. A prototype of the proposed design is fabricated to verify the numerical results of the model. Human motions, such as finger pressing, have been proven to trigger an PEH for driving a batteryless switch [26]. In this study, the PEH is driven by finger pressing and can be applicable to such applications.

This paper is organized as follows. In Section 1, the background of this study is presented. Some selected papers for piezoelectric energy harvesting are reviewed. The working principle of the PEH is mentioned in Section 2. In Section 3, the mathematical model for calculating the dynamic behaviors of the PEH are presented in detail. The fabrication and experimental setup of the PEH prototype are presented in Section 4. The discussions on the numerical and experimental results are presented in Section 5. The conclusions and findings of this study are summarized in Section 6.

#### **2. Working Principle of the PEH**

Figure 1a shows a conceptual drawing of the PEH, which includes a bimorph piezoelectric cantilever beam and three magnets. The magnetizations of the magnets are marked by N and S. The moving directions of the magnet A and the beam are also marked in Figure 1a. Figure 1b shows

the side view and front view of the PEH. The magnet B is mounted at the beam's tip. At its initial position, the beam magnet B is attracted on the moving magnet A. When magnet A moves along *y*-direction, the magnetic force acting on magnet B become smaller. At a critical position *u*0, the elastic force of the beam equals the magnetic force, and the beam magnet B is about to separate from magnet A, as shown in Figure 1c. The upward motion of the beam is driven by an upward force *Fb*, which is the resultant of the magnetic forces from magnets A and C, and the elastic force of the beam bending. It should be noted that at the critical position *u*0, the beam magnetic force could produce a minor torque on the beam. It is assumed that the torque is small so that the induced torsion deformation can be neglected. During its upward-moving, the beam magnet B has a position *uz* and a velocity . *uz*. Meanwhile, the magnet A has a position *uy* and a velocity . *uy*, as illustrated in Figure 1d. Finally, the beam magnet B collides magnet C. After the impact, a transient vibration in the beam occurs and the electric energy is then extracted from the vibration due to the piezoelectric effect.

**Figure 1.** The conceptual drawing of the piezoelectric energy harvester (PEH). (**a**) The iso-view of the PEH at the initial status; (**b**) The side view and front view of the PEH at the initial status; (**c**) The PEH at the time that magnet B separates from magnet A; (**d**) The beam deflection before impact.

In the proposed PEH, the magnet B is moved by finger pressing. It means that the PEH is driven by human motion and could be applicable for driving a batteryless switch [26].

#### **3. Dynamic Model of the PEH**

As shown in Figure 2, the piezoelectric beam is composed of a pair of PZT (lead zirconate titanate) and a middle metal shim. The beam is divided by 7 sub-beams. According to the motion sequence of the PEH, the analysis is divided into two phases. In the first phase, we consider the dynamic responses of the beam motion before the impact. In the second phase, the transient vibration after the impact is solved. In the rest of this section, the mathematical procedures of the two phases will be derived separately.

**Figure 2.** The beam model of the piezoelectric beam.

#### *3.1. Analysis before Impact*

The first phase of the model is to consider the beam motion before impact. In this phase, the acceleration is much smaller than that after impact so that we only consider the fundamental vibration mode of vibration, which is assumed to be the same as the static beam deflection subjected to a concentrated force at its end. Based on the Euler beam theory, the static beam deflection can be easily determined by the flexural formula. In a general form, the deflection can be written as

$$w\_k(\mathbf{x}\_k, t) = w\_k^{(1)}(\mathbf{x}\_k, t) = \frac{w\_b(t)}{l!l\_b} W\_k^{(1)}(\mathbf{x}\_k), \ k = 0, 1, \cdots \text{-6} \tag{1}$$

where the superscript (1) denotes the first phase, *<sup>w</sup>*(1) *<sup>k</sup>* (*xk*, *<sup>t</sup>*) denotes the deflection of the *<sup>k</sup>*th sub-beam before impact, *<sup>W</sup>*(1) *<sup>k</sup>* (*xk*) are the normalized deflection of the *<sup>k</sup>*th sub-beam under a unit force applied at the beam tip, and *Ub* is the normalized deflection at the beam tip, i.e., *Ub* <sup>=</sup> *<sup>W</sup>*(1) <sup>6</sup> (*x*<sup>6</sup> = *L*6). The expressions for *<sup>W</sup>*(1) *<sup>k</sup>* (*xk*) is given in the Appendix A. It is seen that *wb* actually denotes the deflection at the beam tip, i.e., *wb*(*t*) = *<sup>w</sup>*(1) <sup>6</sup> (*x*<sup>6</sup> = *L*6).

By neglecting the small tilt angle θ of beam magnet B and the magnetic force along *x* and *y* directions, the force *Fa*-*<sup>b</sup>* along *z* direction acting on magnet B from magnet A can be written in the form [27]

$$F\_{a-b} = \frac{M\_{\text{fl}}M\_{\text{b}}}{4\pi\mu\_0} (\phi\_1 + \phi\_2 + \phi\_3 + \phi\_4) \tag{2}$$

where *Ma* and *Mb* denote the magnetizations of the magnets A and B, respectively, μ0 = <sup>4</sup><sup>π</sup> <sup>×</sup> <sup>10</sup>−<sup>7</sup> N/A2 denotes the permeability of the air, and

$$\phi\_1 = -\sum\_{i=0}^{1} \sum\_{j=0}^{1} \sum\_{k=0}^{1} \sum\_{l=0}^{1} \sum\_{p=0}^{1} \sum\_{q=0}^{1} \left[ \mu\_{ij} w\_{pq} \ln \left( \sqrt{u\_{ij}^2 + v\_{kl}^2 + w\_{pq}^2} - u\_{ij} \right) \right] \tag{3}$$

$$\phi\_2 = -\sum\_{i=0}^{1} \sum\_{j=0}^{1} \sum\_{k=0}^{1} \sum\_{l=0}^{1} \sum\_{p=0}^{1} \sum\_{q=0}^{1} \left[ \upsilon\_{kl} w\_{pq} \ln \left( \sqrt{u\_{ij}^2 + v\_{kl}^2 + w\_{pq}^2} - \upsilon\_{kl} \right) \right] \tag{4}$$

$$\phi\_3 = \sum\_{i=0}^{1} \sum\_{j=0}^{1} \sum\_{k=0}^{1} \sum\_{l=0}^{1} \sum\_{p=0}^{1} \sum\_{q=0}^{1} \left| u\_{ij} v\_{kl} \tan^{-1} \left( \frac{u\_{ij} v\_{kl}}{w\_{pq} \sqrt{u\_{ij}^2 + v\_{kl}^2 + w\_{pq}^2}} \right) \right| \tag{5}$$

$$\phi\_4 = -\sum\_{i=0}^{1} \sum\_{j=0}^{1} \sum\_{k=0}^{1} \sum\_{l=0}^{1} \sum\_{p=0}^{1} \sum\_{q=0}^{1} \left[ w\_{pq} \sqrt{u\_{ij}^2 + v\_{kl}^2 + w\_{pq}^2} \right] \tag{6}$$

In Equations (3)–(6), the parameters *uij*, *vkl*, and *wpq* are given by

$$\iota\_{i\dot{j}} = \alpha + (-1)^{l} L\_{\mathfrak{a}} - (-1)^{l} L\_{\mathfrak{b}\prime} \iota\_{\mathfrak{t}\prime} \dot{\mathfrak{z}} = 0, \ 1 \tag{7}$$

$$
\omega\_{kl} = \beta + (-1)^l \mathcal{W}\_{\mathfrak{d}} - (-1)^k \mathcal{W}\_{\mathfrak{h}} \text{ } k \text{ } l = \text{ 0, 1} \tag{8}
$$

$$w\_{pq} = \gamma + (-1)^q H\_a - (-1)^p H\_{b\prime} \ p \; \; q \; = \; 0 \; \; 1 \tag{9}$$

where (*La*, *Ha*, *Wa*) and (*Lb*, *Hb*, *Wb*) be the dimensions of the two magnets, (α, β, γ) is relative position vector from the center of magnet A to the center of magnet B. As shown in Figure 1d, it is seen that the components of the position vector are α = 0, β = −*uy* and γ = *uz* in this particular case.

Following a very similar procedure, *Fc*-*b*, the *z*-component of magnetic force acting on the beam magnet B from magnet C, can be also calculated. The resultant force of the two magnetic forces is *Fb* = *Fa*-*<sup>b</sup>* + *Fc*-*b*. Initially, the magnetic force *Fb* is greater than the elastic force of the beam and the magnet B is attracted by and at rest on the magnet A. When the magnet A begins to move upward along *y* direction, the magnetic force *Fb* decreases. At a critical status *uy* = *u*0, the beam begins to move upward. During the upward motion, the kinetic energy *<sup>T</sup>*(1) elastic internal energy *<sup>V</sup>*(1) *<sup>e</sup>* and electrostatic internal energy *<sup>V</sup>*(1) *<sup>s</sup>* can be written in the following forms:

$$T^{(1)} = \frac{1}{2} \sum\_{k=0}^{6} \int\_{0}^{L\_k} (\rho A)\_k \left(\dot{w}\_k^{(1)}\right)^2 d\mathbf{x}\_k \tag{10}$$

$$V\_c^{(1)} = \frac{1}{2} \sum\_{k=0}^{5} \int\_0^{l\_k} (EI)\_k \left(\frac{d^2 w\_k^{(1)}}{dx\_k^2}\right)^2 dx\_k - \frac{bd\_{31} (t\_m + t\_p) V}{4s\_{11}^E} \sum\_{k=1}^4 \int\_0^{l\_k} \frac{\partial^2 w\_k^{(1)}}{\partial x^2} dx\_k \tag{11}$$

$$V\_s^{(1)} = \frac{bd\_{31} \binom{t\_m t\_p}{} V^{(1)} \frac{4}{\
u}}{4s\_{11}} \sum\_{k=1}^4 \int\_0^{L\_k} \frac{\partial^2 w\_k^{(1)}}{\partial \mathbf{x}\_k^2} d\mathbf{x}\_k + \left(\frac{\varepsilon\_{33} s\_{11} - d\_{31}^2}{s\_{11}}\right) \frac{b (L\_1 + L\_2 + L\_3 + L\_4) \left(V^{(1)}\right)^2}{4t\_p} \tag{12}$$

where *V*(1) denotes the output voltage of the PEH before impact, (*EI*)*<sup>k</sup>* is the bending rigidity of the *k*th sub-beam, *b* is the width of the beam, *d*<sup>31</sup> is the piezoelectric constant, *s*<sup>11</sup> is the compliance of the piezoelectric, ε<sup>33</sup> is the dielectric constant of the piezoelectric, *tm* is the thickness of the metal shim, and *tp* is the thickness of the piezoelectric. The Lagrange equation for the 1-D motion can be written as

$$\frac{d}{dt}\left(\frac{\partial L^{(1)}}{\partial \dot{w}\_b}\right) - \frac{\partial L^{(1)}}{\partial w\_b} = F\_b \tag{13}$$

where the Lagragian *<sup>L</sup>*(1) is defined by *<sup>L</sup>*(1) <sup>=</sup> *<sup>T</sup>*(1) <sup>−</sup> *<sup>V</sup>*(1) *<sup>e</sup>* <sup>+</sup> *<sup>V</sup>*(1) *<sup>s</sup>* . The current generated by the PEH is

$$\dot{q}^{(1)} = -\frac{bd\_{31}(t\_m + 2t\_p)}{2s\_{11}} \sum\_{k=1}^{4} \int\_0^{L\_k} \frac{\partial^2 \dot{w}\_k^{(1)}}{\partial x\_k^2} d\mathbf{x}\_k + \frac{b\left(d\_{31}^2 - s\_{11}\varepsilon\_{33}\right)(L\_1 + L\_2 + L\_3 + L\_4)\dot{V}^{(1)}}{2t\_pS\_{11}} \tag{14}$$

If the PEH is connected to an external resistance *R*, the equation of circuit can be obtained by the Ohm's law, i.e.,

$$i^{(1)} = \frac{V^{(1)}}{R} \tag{15}$$

Substituting Equations (10)–(12) into Equation (13) and substituting Equation (14) into Equation (15), the equation of motion and equation of circuit before impact are

$$
\dot{m}\_{\alpha\eta}\ddot{w}\_b + k\_{\alpha\eta}w\_b + \alpha V = F(u\_{y'}, u\_z) \tag{16}
$$

*Sensors* **2020**, *20*, 6170

$$
\eta \alpha \dot{w}\_b - \mathbb{C}\_p \dot{V} = \frac{V}{R} \tag{17}
$$

where

$$u\_z(t) = w\_6^{(s)}(\mathbf{x}\_6 = L\_6/2, \mathbf{t}) + \mathbf{h} + \frac{H\_d}{2} \tag{18}$$

$$\alpha = -\frac{bd\_{31}\Big(t\_m t\_p\Big)}{4s\_{11} l I\_b} \sum\_{k=1}^4 \int\_0^{L\_k} \frac{\partial^2 \mathcal{W}\_k^{(s)}(\mathbf{x}\_k)}{\partial \mathbf{x}\_k^2} d\mathbf{x}\_k \tag{19}$$

$$
\eta = \frac{2t\_m + 4t\_p}{t\_m t\_p} \tag{20}
$$

$$\mathcal{C}\_{p} = \frac{b(L\_1 + L\_2 + L\_3 + L\_4)}{2t\_p} \Big(\varepsilon\_{33} - \frac{d\_{31}^2}{s\_{11}}\Big) \tag{21}$$

Assume that the magnet A is moving at a constant velocity *va*. Then the position of the magnet A is

$$
\mu\_y = \mu\_0 + \upsilon\_a(t - t\_0) \tag{22}
$$

where *t*<sup>0</sup> is the time that the beam magnet B is about to separate from the magnet A. The initial condition at *t* = *t*<sup>0</sup> is

$$\dot{u}\_y = u\_{0,}\ w\_b = -h - \frac{H\_a}{2},\ \dot{w}\_b = 0\tag{23}$$

By solving Equations (16) and (17) with the initial conditions in Equation (23), the beam tip deflection *wb* and output voltage *V*(1) before impact can be solved. By substituting the solved *wb* into Equation (1), the deflection curve of the piezoelectric beam before impact can be determined.

#### *3.2. Analysis after Impact*

At the instance of the impact, the deflection at the beam tip is *wb* = *h* − *Hb*/2. By the use of Equation (1), the deflection of the beam at this instance is

$$w\_k^{(1)}(\mathbf{x}\_k, t=0) = \frac{2h - H\_b}{2Ll\_b} \mathcal{W}\_k^{(1)}(\mathbf{x}\_k), \ k = 0, 1, \cdots \cdot 6\tag{24}$$

Figure 3a shows the illustration of the deflection curve at the impact and after the impact. By introduce the dynamic displacement *<sup>w</sup>*(2) *<sup>k</sup>* (*xk*,*t*) after the impact, the deflection of the beam *<sup>w</sup>*(2) *k* after the impact can be written as

$$w\_k(\mathbf{x}\_{k'}t) = w\_k^{(2)}(\mathbf{x}\_{k'}t) + w\_k^{(1)}(\mathbf{x}\_{k'}t = 0), \; k = 0, 1, \cdots \text{6} \tag{25}$$

Note that the second term in the right-hand-side of Equation (25) is independent of time so that the velocity is . *wk*(*xk*,*t*) = . *<sup>w</sup>*(2) *<sup>k</sup>* (*xk*,*t*). The dynamic displacement *<sup>w</sup>*(2) *<sup>k</sup>* (*xk*,*t*) can be written in linear combination of the interpolation functions *Nk*1(*xk*) to *Nk*4(*xk*), *k* = 0, 1, ... , 6, i.e.,

$$w\_k^{(2)}(\mathbf{x}\_k, t) = q\_k^{(2)}(t)\mathcal{N}\_{k1}(\mathbf{x}\_k) + \theta\_k^{(2)}(t)\mathcal{N}\_{k2}(\mathbf{x}\_k) + q\_{k+1}^{(2)}(t)\mathcal{N}\_{k3}(\mathbf{x}\_k) + \theta\_{k+1}^{(2)}(t)\mathcal{N}\_{k4}(\mathbf{x}\_k), \ k = 0, 1, \cdots, 6 \tag{26}$$

where *q* (2) *<sup>k</sup>* (*t*) and <sup>θ</sup>(2) *<sup>k</sup>* (*t*) denote the dynamic displacement and dynamic rotation at the *k*th node for the dynamic term. The expressions for the interpolation functions are given in the Appendix A.

The magnet B may rebound after it impacts the magnet C. To estimate the rebound displacement, consider a simplified model that the magnet B approaches magnet C and collision occurs between them. In this simplified model, the beam is neglected. The rebound velocity of magnet B after impact is *e* . *wb wb*=*h*−*Hb*/2, where . *wb wb*=*h*−*Hb*/2 denotes the impact velocity of magnet B at the instance just before the impact and *e* is the coefficient of restitution during the impact. By knowing the mass of

the magnet, the kinetic energy of the magnet B can be calculated according to the impact velocity. During the rebound, a negative work is done by the magnetic force. By using the work–kinetic energy principle, the rebound displacement can be estimated. It should be noted that the above simple model overestimates the rebound displacement because we neglect the effects from the beam, which has an upward momentum during the impact and reduce the downward rebound displacement. For the case of small rebound displacement, one can assume that the contact point (corner point 1 shown in Figure 3b) remains no separation after impact. The other corner point 2 shown in Figure 3b is allowed to have motion. However, the degree of freedom of this point is constrained by the magnetic force between two magnets. In this model, an equivalent spring *ks* between the magnets B and C is introduced to model this constraint, as shown in Figure 3b. For simplification without loss of generality, we set *t* = 0 at the instance of the impact. In the following, the analysis procedures of the beam after the impact will be derived in detail.

**Figure 3.** (**a**) Schematic of the beam deflection after impact; (**b**) The spring model of the contact.

In the equivalent spring model, the spring force *Fs*(*t*) is given by

$$F\_s(t) = k\_s \delta(t) \tag{27}$$

where *ks* denotes the equivalent spring constant and

$$\delta(t) = h - w\_{\boldsymbol{\theta}}^{(2)}(\mathbf{x}\_{\boldsymbol{\theta}} = 0, t) - w\_{\boldsymbol{\theta}}^{(1)}(\mathbf{x}\_{\boldsymbol{\theta}} = 0, t = 0) \tag{28}$$

As an estimation, *ks* is written in the form

$$k\_s = \frac{dF}{d\delta} \tag{29}$$

To derive the equation of motion, consider the elastic internal energy *<sup>V</sup>*(2) *<sup>e</sup>* , electrostatic internal energy *<sup>V</sup>*(2) *<sup>s</sup>* and kinetic energy. The equations of motion can be derived by Lagrange mechanics:

$$\frac{d}{dt}\left(\frac{\partial L^{(2)}}{\partial \dot{q}\_k^{(2)}}\right) - \frac{\partial L^{(2)}}{\partial q\_k^{(2)}} = 0, \; k = 1, 2, \cdots, 6 \tag{30}$$

$$\frac{d}{dt}\left(\frac{\partial L^{(2)}}{\partial \dot{\theta}\_k^{(2)}}\right) - \frac{\partial L^{(2)}}{\partial \theta\_k^{(2)}} = 0, \; k = 1, 2, \cdots, 5 \tag{31}$$

where *<sup>L</sup>*(2) <sup>=</sup> *<sup>T</sup>*(2) <sup>−</sup> *<sup>V</sup>*(2) *<sup>e</sup>* <sup>−</sup> *<sup>V</sup>*(2) *<sup>e</sup>* . The kinetic energy and electrostatic internal energy for the second phase have the same forms as those mentioned in Equations (10)–(12) by changing superscript (1) by (2). The elastic internal energy after impact is

$$V\_{\varepsilon}^{(2)} = \frac{1}{2} \sum\_{k=0}^{5} (EI)\_k \int\_0^{l\_k} \left(\frac{\partial^2 w\_k^{(2)}}{\partial \mathbf{x}^2}\right)^2 d\mathbf{x}\_k - \frac{bd\_{31} (t\_m + t\_p) V}{2s\_{11}} \sum\_{k=1}^{4} \int\_0^{l\_k} \frac{\partial^2 w\_k^{(2)}}{\partial \mathbf{x}^2} d\mathbf{x}\_k + \frac{1}{2} k\_{4\eta} \left(\mathbb{I} - q\_6^{(2)} - q\_6^{(1)}\right)^2 \tag{32}$$

Similarly, the circuit equation after impact can be derived in the same way. The equations of the motion and equation of circuit after impact can be written by

$$\mathbf{M}\ddot{\mathbf{q}}^{(2)} + \mathbf{C}\dot{\mathbf{q}}^{(2)} + \mathbf{K}\mathbf{q}^{(2)} + \mathbf{a}V^{(2)} = 0\tag{33}$$

$$
\eta \mathbf{a}^T \dot{\mathbf{q}}^{(2)} - \mathbf{C}\_P \dot{\mathbf{V}}^{(2)} = \frac{V^{(2)}}{R} \tag{34}
$$

where **M** is the mass matrix considering the masses of the beam and magnet, **C** is the damping matrix, **K** is the stiffness matrix considering both the elastic force of the beam and the equivalent spring, α is the electromechanical coupling matrix containing the piezoelectric constant, *Cp* is the capacitor of the PZT, *V* is the output voltage, *R* is the resistance, and

$$\mathbf{q}^{(2)} = \begin{bmatrix} q\_1^{(2)} & \theta\_1^{(2)} & \cdots & q\_5^{(2)} & \theta\_5^{(2)} & q\_6^{(2)} \end{bmatrix}^T \tag{35}$$

The expressions for the matrices **M**, **K** and α are given in the Appendix A. The proportional damping model is used in the present analysis, i.e.,

$$\mathbf{C} = \gamma\_1 \mathbf{M} + \gamma\_2 \mathbf{K} \tag{36}$$

where γ<sup>1</sup> and γ<sup>2</sup> are constants and can be determined by solving the following two equations:

$$\frac{\mathcal{V}\_1}{\omega\_1} + \omega\_1 \mathcal{V}\_2 = 2\zeta\_1 \tag{37}$$

$$a\_2 \frac{\mathcal{V}\_1}{\omega\_2} + a\_2 \mathcal{V}\_2 = 2\zeta\_2 \tag{38}$$

where ω*k*, ζ*<sup>k</sup>* (*k* = 1, 2) are the natural frequency and damping ratio, respectively, of the first two vibration modes.

The boundary conditions for the dynamic displacement at *x*<sup>0</sup> = 0 and *x*<sup>6</sup> = *L*<sup>6</sup> after the impact are

$$
\theta\_0^{(2)} = 0, \; \theta\_0^{(2)} = 0, \; \theta\_7^{(2)} = 0 \tag{39}
$$

At the instance that the impact occurs (*t* = 0), the initial conditions for the second phase are

$$\begin{aligned} \left.q\_k^{(2)}(t=0)=0, \left.\theta\_k^{(2)}(t=0)=0, \left.\dot{q}\_k^{(2)}(t=0)\right.\right. &= \left.\frac{1}{L\_b}W\_k^{(1)}\right|\dot{w}\_b\Big|\_{w\_b=h-\frac{H\_b}{2}}\\ \left.\dot{\theta}\_k^{(2)}(t=0)=\frac{1}{L\_b}\frac{d\left|\mathcal{W}\_k^{(2)}\right|}{dx\_b}\Big|\dot{w}\_b\Big|\_{w\_b=h-\frac{H\_b}{2}}\right), \left.V^{(2)}(t=0)=V^{(1)}\right|\_{w\_b=h-\frac{H\_b}{2}} \end{aligned} \tag{40}$$

where . *wb wb*=*h*<sup>−</sup> *Hb* 2 denotes the velocity of the beam tip at the time just before the impact and can be determined from the solution of Equation (13) mentioned in phase 1.

By solving Equations (33) and (34) in conjunction with the boundary conditions Equation (39) and initial conditions Equation (40), the transient responses of displacements and voltage of the beam can be obtained.

The energy *E*PEH harvested by the PEH can be calculated by

$$E\_{\rm PEH} = \int\_{t\_0}^{0} \frac{V^2}{R} dt + \int\_0^{\Delta t} \frac{V^2}{R} dt \tag{41}$$

where Δ*t* is the time interval for calculating the energy.

#### **4. Fabrication of Prototype and Experimental Setup**

Figure 4a illustrates the PEH prototype. The piezoelectric beam is clamped on a base. The magnets A and C are respectively mounted on two sliding bars, which can slide under the guidance of a pair of guiders. In the initial status, the beam magnet B is attracted by magnet A and is at rest. The upper sliding bar is placed at a position that the magnet C is at the position aligned to magnets A and B. By pushing the sliding bar to move right, the magnet A separates magnet B, and the beam moves upward and finally stopped by the magnet C. Then an impact occurs. In Figure 4b, the two bars are hidden for a better view for piezoelectric beam and magnets. The base, guiders and bars are made by acrylic to avoid interfering the magnetic fields.

**Figure 4.** (**a**) The design of the PEH; (**b**) The design of the PEH with invisible slider bars; (**c**) Piezoelectric beam; (**d**) Prototype of the PEH.

The fabrication and assembly of the piezoelectric beam were provided by Eleceram Technology Co., Taiwan, as shown in Figure 4c. The geometric parameters of the beam and magnets are listed in Table 1. Two identical NdFeB magnets with a size of 10 <sup>×</sup> 5 <sup>×</sup> 0.8 mm<sup>3</sup> were adhered at the beam tip to form the magnet B, which has an equivalent size of 10 <sup>×</sup> 5 <sup>×</sup> 1.6 mm3. Another two identical NdFeB magnets serve as magnets A and C, each of which has a size of 10 <sup>×</sup> <sup>10</sup> <sup>×</sup> 10 mm3. The geometric parameters of the magnets are also summarized in Table 1. Figure 4d shows the assembled prototype.


**Table 1.** The geometric parameters of the beam and magnets.

Figure 5 shows the setup for measuring the magnetic forces from one magnet to another for various distances between two magnets. In the experiment, a magnet was fixed on the fixture. A second magnet was placed on a jig that can freely slide along a guider. A force gauge (DigiTech DTG-10) withstood the second magnet at a gap between the two magnets and measured the magnetic force. Three magnet pairs A-A, B-B and A-B were considered in the magnetic force measurements for various gaps between the two magnets.

**Figure 5.** The experimental setup for the magnetic force measurement.

For the magnet pair A-A experiment, the magnetic force *F* (exp) *<sup>a</sup>*−*a*,*<sup>k</sup>* is measured for *<sup>k</sup>*th gap. Based on the 3-D magnetic force model described in Equation (2), the magnetic force *<sup>F</sup>*(model) *<sup>a</sup>*−*a*,*<sup>k</sup>* of the magnet pair A-A can be written in the form:

$$F\_{a-a,k} = \frac{M\_a^2}{4\pi\mu\_0} (\phi\_{1,k} + \phi\_{2,k} + \phi\_{3,k} + \phi\_{4,k}) \tag{42}$$

where the subscript *k* denotes the parameters calculated according to the *k*th gap. According to the least square method, the magnetic magnetization *Ma* of the magnet A can be determined solving the equation:

$$\frac{d}{dQ\_{a-a}}\sum\_{k=1}^{N} \left( F\_{a-a,k}^{(\text{exp})} - F\_{a-a,k} \right)^2 = 0 \tag{43}$$

where *Qa*−*<sup>a</sup>* = *<sup>M</sup>*<sup>2</sup> *<sup>a</sup>* and *N* denotes the number of gaps for the measurement. Following a similar procedure, the magnetization *Mb* can be also determined from the magnet pair B-B experiment.

To investigate the energy generated by the PEH, an oscilloscope (Tektronix DPO 4054B) was used to monitor the output voltage of the resistor connected to the PEH, as shown in Figure 6a. Because the motion of the sliding bar is moved by hand, the velocities of the magnet A change for different tests. The motion of the sliding bar (or the magnet A) was monitored by a high-speed camera (OLYMPUS i-SPEED 3) with a frame rate of 5000 Hz, as shown in Figure 6b. By analyzing the photos frame by frame, the positions of the sliding bar and the beam tip can be determined. The velocities of the beam magnet B can be calculated accordingly.

**Figure 6.** The experimental setup of the PEH: (**a**) voltage measurement of the PEH coupled with a resistance; (**b**) motion measurement of the PEH by the highspeed camera; (**c**) voltage measurement of the PEH coupled with a diode-bridge, capacitor, and LED.

Figure 6c shows experimental setup for light LED. Four diodes (1N4004) were connected to form a bridge. The AC voltages were rectified to DC signals and the electric charges were stored in a capacitor *Cst* of 10 μF. An LED was connected to the capacitor in parallel. In the experiments, the first impact was driven by magnet A and the collision between the magnets B and C. Put the magnet A back to the initial position and move the magnet C. Then the beam moved to magnet A and the second impact occurred. The process was repeated ten times and produced ten impacts. For each impact, the capacitor was charged and the voltage across the capacitor was boosted. After ten impacts, the switch S was switched to connect the LED loop. The capacitor discharged and a current flowed through and lighted up the LED. The oscilloscope was used to monitor the voltage across the capacitor during the charging and discharging in order to evaluate the energy stored in the capacitor. A multimeter (Fluke 189) was used

to capture the peak current during the discharging in order to evaluate the maximum instantaneous power for lighting the LED.

#### **5. Results and Discussions**

In this section, the experimental and numerical results are presented to demonstrate the performance of the PEH. The material properties used in the dynamic model are listed in Table 2.


**Table 2.** The material properties used in the model.

#### *5.1. Measurements of Magnetic Forces*

The variation of the magnetic forces for various gaps are shown in Figure 7. Note that the magnetizations *Ma*, *Mb*, and *Mc* for the magnets A, B, and C are respectively unknown. By using the least square method described in Equation (43), the magnetizations *Ma* and *Mb* are determined to be 0.9537 T and 0.5147 T, respectively. In Figure 7, it is seen that the regression curves for magnet pairs A-A and B-B agree well with the measurement data. These two magnetizations are used to calculate the magnetic forces for the magnet pair A-B. In Figure 7, it is seen that the computation results for magnetic force between magnets A and B are validated by experiments. It indicates that Equation (2) can used to predict the magnetic force at high accuracy.

**Figure 7.** The comparisons of magnetic forces by measurements and model for magnet pairs A-A, B-B, and A-B for various gaps.

#### *5.2. E*ff*ects of the Velocity of the Magnet A*

The PEH was driven by finger pressing. The velocity *va* of the magnet A varies for different tests. To investigate the effects of the velocity *va* on the motion of the beam and output voltage, five tests to trigger the PEH have been tested. Figure 8 shows the measured beam tip deflections *wb* for the first test. The deflection results show that the beam tip accelerated during its upward motion before impact. The measured positions *uy* of the magnet A for the first test are also showed in Figure 8. It is seen that the position *uy* of the magnet A exhibits a linear trend over time and the constant velocity assumption described in Equation (22) can be acceptable.

**Figure 8.** The measured deflection *wb* of the beam tip and the position *uy* of the magnet A before impact for the first test.

Table 3 lists the measured relative position *u*<sup>0</sup> between magnets A and B at the instance that the beam tip begins to move upward, and the measured velocities *va* of the magnet A for the five tests. The variation of relation position *u*<sup>0</sup> for different tests is small. Meanwhile, the velocity *va* varies from different tests because the motion of the magnet A was driven by hand. Among the five tests, the maximum, minimum and average values of *va* are 0.1337 m/s, 0.07763 m/s, and 0.09731 m/s, respectively. The measured impact velocity, which is defined as the beam tip velocity . *wb wb*=*h*<sup>−</sup> *Hb* 2 at the time just before impact, for the five tests are also listed in Table 3. It is observed that the relation between impact velocity and the moving magnet velocity *va* is insignificant. The average impact velocity for the five tests is 5.106 m/s. Assume that the coefficient of restitution during the impact is *e* = 0.5, then the rebound velocity of the magnet B is 2.553 m/s. By knowing the mass of 0.584 g for the magnet B, the kinetic energy after impact is 1.955 mJ. According to the magnetic force measurements, the magnetic attraction force is approximately 3 N when the magnet B is close to magnet C. By using the energy balance, the rebound displacement is 0.65 mm, which is quite small by comparing the whole moving distance. As mentioned in Section 3.2, the calculation of the rebound displacement is overestimated. It can be concluded that the small rebound assumption is acceptable. In Figure 8, the measurement data also show the beam tip deflection after impact. It is observed that the beam tip becomes stationary after impact, i.e., no rebound can be observed. Therefore, the model shown in Figure 3b can be acceptable.

Figure 9 shows the time history of the measured beam tip deflection for the five tests. Each test reveals almost identical beam tip deflection curve. It indicates that the beam tip motion is independent of the driving velocity *va*. The model calculation results agree well with the measured data, indicating that the model for the first phase can predict the motion of the beam tip at high accuracy.

Figure 10 shows the time history of magnetic forces applied on magnet B from magnets A and C during the phase before impact. As expected, the magnetic force from magnet C increases by time, while the magnetic force from magnet A is relatively small, especially for the time approaching impact. At the time of 1 ms before impact, the magnetic force from magnet C dominates the beam tip motion. The small force magnetic force from magnet A can be used to explain the insignificant contribution of *va*.

Figure 11 shows the open-circuit voltage for the five tests. The transient responses of the voltages before and after impact keep the same for various *va*. The voltage responses computed from the model are also shown in Figure 11. According to the voltage responses after impact shown in Figure 11, the motion mainly exhibits the fundamental mode although some minor higher order vibration modes are observed. The domination of the fundamental mode can be also seen in model calculations. It is seen that the model agrees with the measurements. By picking the peak and valley points of the transient voltage response, the natural frequency of the vibration can be determined. According to this method, the natural frequencies for the five tests in Figure 11 are 407.9, 403.3, 408.9, 413.6, and 410.9 Hz. The average natural frequency for the five data is 408.9 Hz. The same method can be also applied in the time response obtained by the model and the result is 373.0 Hz.

**Table 3.** Measurement results of *u*0, *va*, and impact velocity before impact for the five tests.


**Figure 9.** The comparisons of the deflection *wb* of the beam tip before impact obtained by experiments and model.

**Figure 10.** The computed magnetic forces by 3-D magnetic force model before impact.

**Figure 11.** The comparisons of the open-circuit voltage responses of the PEH obtained by experiments and the dynamic model.

#### *5.3. Energy Measurements of the PEH*

Figure 12 shows the measured maximum *V*max and energy *E*PEH harvested by the PEH for various external resistances. For small external resistance, both *V*max and *E*PEH are small. The energy increases with the increase of *R* when *R* is less than 20 *k*Ω. A maximum energy of 0.4045 mJ occurs at the optimum resistance *R*opt = 15 *k*Ω. For a large *R*, *V*max reaches a stable value and the average power becomes small. In Figure 12, the measured voltage and energy are also compared with the results computed by the model. The comparisons show that the results of the model agree well with the experimental data.

**Figure 12.** (**a**) The maximum output voltage of the PEH for various external resistance *R*; (**b**) The energy harvested by the PEH for various external resistance *R*.

For the lighting up LED experiment, two tests were conducted and the voltages across the capacitor, as shown in Figure 6c, are plotted in Figure 13. For each test, the voltage was boosted to 6.2 V after ten impacts and the energy stored in the capacitor was 19.22 μJ. When the circuit was switched, the capacitor discharged then a current flowed through and lighted up the LED. The voltage drops during the discharging for the two tests were 3.44 V and 3.52 V. The maximum currents during the discharging for the two tests were 450 μA and 476 μA. The maximum instantaneous powers during the discharging for the two tests were 0.774 mW and 0.838 mW. The experiments showed that the energy stored in the capacitor can light up the LED.

**Figure 13.** The measured voltage of the capacitor in the lighting up LED experiment.

#### **6. Conclusions**

In this study, the impact-driven PEH in the magnetic field is presented. The multi-DOF mathematical model was developed for solving the dynamic and vibration behaviors of the piezoelectric beam under the effects of magnetic fields. The prototype was also fabricated and the performance of the PEH was measured in detail. The conclusions of this study are summarized as follows:


**Author Contributions:** Conceptualization, C.-D.C., Y.-H.W., and P.-W.S.; formal analysis, C.-D.C. and Y.-H.W.; methodology Y.-H.W. and P.-W.S.; validation, C.-D.C.; writing—original draft, C.-D.C. and Y.-H.W.; writing—review and editing, C.-D.C. All authors have read and agreed to the published version of the manuscript.

**Funding:** The authors are grateful to the Ministry of Science and Technology, Taiwan, for the financial support through grant MOST 107-2221-E-006-154-MY2.

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

#### **Appendix A**

The normalized deflection *<sup>W</sup>*(1) *<sup>k</sup>* (*xk*) appeared in Equation (1) is given by

$$\mathcal{W}\_{k}^{(1)}(\mathbf{x}\_{k}) = \begin{cases} \frac{1}{2\left(\delta I\right)\_{0}} \Big( \frac{1}{2} L\_{6} + \sum\_{l=0}^{5} L\_{l} \Big) \mathbf{x}\_{0}^{2} - \frac{1}{6\left(\delta I\right)\_{0}} \mathbf{x}\_{0}^{3}, \text{ for } k = 0\\\frac{1}{2\left(\delta I\right)\_{k}} \Big( \frac{1}{2} L\_{6} + \sum\_{l=k}^{5} L\_{l} \Big) \mathbf{x}\_{k}^{2} - \frac{1}{6\left(\delta I\right)\_{k}} \mathbf{x}\_{k}^{3} + \left. \left( \frac{d\mathcal{W}\_{k-1}^{(1)}}{\mathrm{d}\mathbf{x}\_{k-1}} \right|\_{\mathbf{x}\_{k-1} = \mathbf{I}\_{k-1}} \right) \mathbf{x}\_{k} + \left. \mathcal{W}\_{k-1}^{(1)} \right|\_{\mathbf{x}\_{k-1} = \mathbf{I}\_{k-1}}, \text{ for } k = 1, 2, \cdots, 5\\\mathcal{W}\_{6}^{(1)}(\mathbf{x}\_{6}) = \left( \left. \frac{d\mathcal{W}\_{5}^{(1)}}{\mathrm{d}\mathbf{x}\_{5}} \right|\_{\mathbf{x}\_{5} = \mathbf{I}\_{5}} \right) \mathbf{x}\_{6} + \left. \mathcal{W}\_{5}^{(1)} \right|\_{\mathbf{x}\_{5} = \mathbf{I}\_{5}}, \text{ for } k = 6 \end{cases}$$

The interpolation functions appeared in Equation (26) is given by

$$N\_{k1} = 1 - \frac{3\mathbf{x}\_k^2}{L\_k^2} + \frac{2\mathbf{x}\_k^3}{L\_k^3}, \\ N\_{k2} = \mathbf{x}\_k - \frac{2\mathbf{x}\_k^2}{L\_k} + \frac{\mathbf{x}\_k^3}{L\_k^2}, \\ N\_{k3} = \frac{3\mathbf{x}\_k^2}{L\_k^2} - \frac{2\mathbf{x}\_k^3}{L\_k^3}, \\ N\_{k4} = -\frac{\mathbf{x}\_k^2}{L\_k} + \frac{\mathbf{x}\_k^3}{L\_k^2}, \\ k = 0, 2, \dots, 6$$

The element *Mi*-*<sup>j</sup>* (*i*, *j* = 1, 2, ... , 11) at *i*th row and *j*th column of the mass matrix **M** appeared in Equation (33) is given by

*<sup>M</sup>*1−<sup>1</sup> <sup>=</sup> <sup>13</sup> 35 *L*0(ρ*A*)<sup>0</sup> + *L*1(ρ*A*)<sup>1</sup> , *<sup>M</sup>*1−<sup>2</sup> <sup>=</sup> *<sup>M</sup>*2−<sup>1</sup> <sup>=</sup> <sup>−</sup> <sup>11</sup> <sup>210</sup> *L*2 <sup>0</sup>(ρ*A*)<sup>0</sup> <sup>−</sup> *<sup>L</sup>*<sup>2</sup> <sup>1</sup>(ρ*A*)<sup>1</sup> , *<sup>M</sup>*1−<sup>3</sup> <sup>=</sup> *<sup>M</sup>*3−<sup>1</sup> <sup>=</sup> <sup>9</sup> <sup>70</sup>*L*1(ρ*A*)1, *<sup>M</sup>*1−<sup>4</sup> <sup>=</sup> *<sup>M</sup>*4−<sup>1</sup> <sup>=</sup> <sup>−</sup> <sup>13</sup> 420*L*<sup>2</sup> <sup>1</sup>(ρ*A*)1, *<sup>M</sup>*2−<sup>2</sup> <sup>=</sup> <sup>1</sup> <sup>105</sup> *L*3 <sup>0</sup>(ρ*A*)<sup>0</sup> <sup>+</sup> *<sup>L</sup>*<sup>3</sup> <sup>1</sup>(ρ*A*)<sup>1</sup> , *<sup>M</sup>*2−<sup>3</sup> <sup>=</sup> *<sup>M</sup>*3−<sup>2</sup> <sup>=</sup> <sup>13</sup> 420*L*<sup>2</sup> <sup>1</sup>(ρ*A*)1, *<sup>M</sup>*2−<sup>4</sup> <sup>=</sup> *<sup>M</sup>*4−<sup>2</sup> <sup>=</sup> <sup>−</sup> <sup>1</sup> 140*L*<sup>3</sup> <sup>1</sup>(ρ*A*)1, *<sup>M</sup>*3−<sup>3</sup> <sup>=</sup> <sup>13</sup> 35 *L*1(ρ*A*)<sup>1</sup> + *L*2(ρ*A*)<sup>2</sup> , *<sup>M</sup>*3−<sup>4</sup> <sup>=</sup> *<sup>M</sup>*4−<sup>3</sup> <sup>=</sup> <sup>−</sup> <sup>11</sup> <sup>210</sup> *L*2 <sup>1</sup>(ρ*A*)<sup>1</sup> <sup>−</sup> *<sup>L</sup>*<sup>2</sup> <sup>2</sup>(ρ*A*)<sup>2</sup> , *<sup>M</sup>*3−<sup>5</sup> <sup>=</sup> *<sup>M</sup>*5−<sup>3</sup> <sup>=</sup> <sup>9</sup> <sup>70</sup>*L*2(ρ*A*)2, *<sup>M</sup>*3−<sup>6</sup> <sup>=</sup> *<sup>M</sup>*6−<sup>3</sup> <sup>=</sup> <sup>−</sup> <sup>13</sup> 420*L*<sup>3</sup> <sup>2</sup>(ρ*A*)2, *<sup>M</sup>*4−<sup>4</sup> <sup>=</sup> <sup>1</sup> <sup>105</sup> *L*3 <sup>1</sup>(ρ*A*)<sup>1</sup> <sup>+</sup> *<sup>L</sup>*<sup>3</sup> <sup>2</sup>(ρ*A*)<sup>2</sup> , *<sup>M</sup>*4−<sup>5</sup> <sup>=</sup> *<sup>M</sup>*5−<sup>4</sup> <sup>=</sup> <sup>13</sup> 420*L*<sup>2</sup> <sup>2</sup>(ρ*A*)2, *<sup>M</sup>*4−<sup>6</sup> <sup>=</sup> *<sup>M</sup>*6−<sup>4</sup> <sup>=</sup> <sup>−</sup> <sup>1</sup> 140*L*<sup>3</sup> <sup>2</sup>(ρ*A*)2, *<sup>M</sup>*5−<sup>5</sup> <sup>=</sup> <sup>13</sup> 35 *L*2(ρ*A*)<sup>2</sup> + *L*3(ρ*A*)<sup>3</sup> , *<sup>M</sup>*5−<sup>6</sup> <sup>=</sup> *<sup>M</sup>*6−<sup>5</sup> <sup>=</sup> <sup>−</sup> <sup>11</sup> <sup>210</sup> *L*2 <sup>2</sup>(ρ*A*)<sup>2</sup> <sup>−</sup> *<sup>L</sup>*<sup>2</sup> <sup>3</sup>(ρ*A*)<sup>3</sup> , *<sup>M</sup>*5−<sup>7</sup> <sup>=</sup> *<sup>M</sup>*7−<sup>5</sup> <sup>=</sup> <sup>9</sup> <sup>70</sup>*L*3(ρ*A*)3, *<sup>M</sup>*5−<sup>8</sup> <sup>=</sup> *<sup>M</sup>*8−<sup>5</sup> <sup>=</sup> <sup>−</sup> <sup>13</sup> 420*L*<sup>2</sup> <sup>3</sup>(ρ*A*)3, *<sup>M</sup>*6−<sup>6</sup> <sup>=</sup> <sup>1</sup> <sup>105</sup> *L*2 <sup>2</sup>(ρ*A*)<sup>2</sup> <sup>+</sup> *<sup>L</sup>*<sup>2</sup> <sup>3</sup>(ρ*A*)<sup>3</sup> , *<sup>M</sup>*6−<sup>7</sup> <sup>=</sup> *<sup>M</sup>*7−<sup>6</sup> <sup>=</sup> <sup>13</sup> 420*L*<sup>2</sup> <sup>3</sup>(ρ*A*)3, *<sup>M</sup>*6−<sup>8</sup> <sup>=</sup> *<sup>M</sup>*8−<sup>6</sup> <sup>=</sup> <sup>−</sup> <sup>1</sup> 140*L*<sup>2</sup> <sup>3</sup>(ρ*A*)3, *<sup>M</sup>*7−<sup>7</sup> <sup>=</sup> <sup>13</sup> 35 *L*3(ρ*A*)<sup>3</sup> + *L*4(ρ*A*)<sup>4</sup> , *<sup>M</sup>*7−<sup>8</sup> <sup>=</sup> *<sup>M</sup>*8−<sup>7</sup> <sup>=</sup> <sup>−</sup> <sup>11</sup> <sup>210</sup> *L*2 <sup>3</sup>(ρ*A*)<sup>3</sup> <sup>−</sup> *<sup>L</sup>*<sup>2</sup> <sup>4</sup>(ρ*A*)<sup>4</sup> , *<sup>M</sup>*7−<sup>9</sup> <sup>=</sup> *<sup>M</sup>*9−<sup>7</sup> <sup>=</sup> <sup>9</sup> <sup>70</sup>*L*4(ρ*A*)4, *<sup>M</sup>*7−<sup>10</sup> <sup>=</sup> *<sup>M</sup>*10−<sup>7</sup> <sup>=</sup> <sup>−</sup> <sup>13</sup> 420*L*<sup>2</sup> <sup>4</sup>(ρ*A*)4, *<sup>M</sup>*8−<sup>8</sup> <sup>=</sup> <sup>1</sup> <sup>105</sup> *L*3 <sup>3</sup>(ρ*A*)<sup>3</sup> <sup>+</sup> *<sup>L</sup>*<sup>3</sup> <sup>4</sup>(ρ*A*)<sup>4</sup> , *<sup>M</sup>*8−<sup>9</sup> <sup>=</sup> *<sup>M</sup>*9−<sup>8</sup> <sup>=</sup> <sup>13</sup> 420*L*<sup>2</sup> <sup>4</sup>(ρ*A*)4, *<sup>M</sup>*8−<sup>10</sup> <sup>=</sup> *<sup>M</sup>*10−<sup>8</sup> <sup>=</sup> <sup>−</sup> <sup>1</sup> 140*L*<sup>3</sup> <sup>4</sup>(ρ*A*)4, *<sup>M</sup>*9−<sup>9</sup> <sup>=</sup> <sup>13</sup> 35 *L*4(ρ*A*)<sup>4</sup> + *L*5(ρ*A*)<sup>5</sup> , *<sup>M</sup>*9−<sup>10</sup> <sup>=</sup> *<sup>M</sup>*10−<sup>9</sup> <sup>=</sup> <sup>−</sup> <sup>11</sup> <sup>210</sup> *L*2 <sup>4</sup>(ρ*A*)<sup>4</sup> <sup>−</sup> *<sup>L</sup>*<sup>2</sup> <sup>5</sup>(ρ*A*)<sup>5</sup> , *<sup>M</sup>*9−<sup>11</sup> <sup>=</sup> *<sup>M</sup>*11−<sup>9</sup> <sup>=</sup> *<sup>L</sup>*5(13*L*5+54*L*6)(ρ*A*)<sup>5</sup> <sup>420</sup>*L*<sup>6</sup> , *<sup>M</sup>*10−<sup>10</sup> <sup>=</sup> <sup>1</sup> <sup>105</sup> *L*3 <sup>4</sup>(ρ*A*)<sup>4</sup> <sup>+</sup> *<sup>L</sup>*<sup>3</sup> <sup>5</sup>(ρ*A*)<sup>5</sup> , *<sup>M</sup>*10−<sup>11</sup> <sup>=</sup> *<sup>M</sup>*11−<sup>10</sup> <sup>=</sup> *<sup>L</sup>*<sup>2</sup> <sup>5</sup>(3*L*5+13*L*6)(ρ*A*)<sup>5</sup> <sup>420</sup>*L*<sup>6</sup> , *<sup>M</sup>*11−<sup>11</sup> <sup>=</sup> <sup>1</sup> 105*L*<sup>2</sup> 6 *L*<sup>3</sup> <sup>5</sup> <sup>+</sup> <sup>11</sup>*L*<sup>2</sup> <sup>5</sup>*L*<sup>6</sup> <sup>+</sup> <sup>39</sup>*L*5*L*<sup>2</sup> 6 (ρ*A*)<sup>5</sup> + <sup>35</sup>*L*<sup>3</sup> <sup>6</sup>(ρ*A*)<sup>6</sup> 

The element *Ki*-*<sup>j</sup>* (*i*, *j* = 1, 2, ... , 11) at *i*th row and *j*th column of the stiffness matrix **K** appeared in Equation (33) is given by

It is noted that the elements not mentioned above for the mass matrix and stiffness matrix are zero.

,

*<sup>K</sup>*1−<sup>1</sup> <sup>=</sup> <sup>12</sup> (*EI*)<sup>0</sup> *L*3 0 + (*EI*)<sup>1</sup> *L*3 1 , *K*1−<sup>2</sup> = *K*2−<sup>1</sup> = 6 <sup>−</sup>(*EI*)<sup>0</sup> *L*2 0 + (*EI*)<sup>1</sup> *L*2 1 , *<sup>K</sup>*1−<sup>3</sup> <sup>=</sup> *<sup>K</sup>*3−<sup>1</sup> <sup>=</sup> <sup>−</sup>12(*EI*)<sup>1</sup> *L*3 1 *<sup>K</sup>*1−<sup>4</sup> <sup>=</sup> *<sup>K</sup>*4−<sup>1</sup> <sup>=</sup> <sup>6</sup>(*EI*)<sup>1</sup> *L*2 1 , *K*2−<sup>2</sup> = 4 (*EI*)<sup>0</sup> *<sup>L</sup>*<sup>0</sup> <sup>+</sup> (*EI*)<sup>1</sup> *L*1 , *<sup>K</sup>*2−<sup>3</sup> <sup>=</sup> *<sup>K</sup>*3−<sup>2</sup> <sup>=</sup> <sup>−</sup>6(*EI*)<sup>1</sup> *L*2 1 , *<sup>K</sup>*2−<sup>4</sup> <sup>=</sup> *<sup>K</sup>*4−<sup>2</sup> <sup>=</sup> <sup>2</sup>(*EI*)<sup>1</sup> *<sup>L</sup>*<sup>1</sup> , *<sup>K</sup>*3−<sup>3</sup> <sup>=</sup> <sup>12</sup> (*EI*)<sup>1</sup> *L*3 1 + (*EI*)<sup>2</sup> *L*3 2 , *K*3−<sup>4</sup> = *K*4−<sup>3</sup> = 6 <sup>−</sup>(*EI*)<sup>1</sup> *L*2 1 + (*EI*)<sup>2</sup> *L*2 2 , *<sup>K</sup>*3−<sup>5</sup> <sup>=</sup> *<sup>K</sup>*5−<sup>3</sup> <sup>=</sup> <sup>−</sup>12(*EI*)<sup>2</sup> *L*3 2 , *<sup>K</sup>*3−<sup>6</sup> <sup>=</sup> *<sup>K</sup>*6−<sup>3</sup> <sup>=</sup> <sup>6</sup>(*EI*)<sup>2</sup> *L*2 2 , *K*4−<sup>4</sup> = 4 (*EI*)<sup>1</sup> *<sup>L</sup>*<sup>1</sup> <sup>+</sup> (*EI*)<sup>2</sup> *L*2 , *<sup>K</sup>*4−<sup>5</sup> <sup>=</sup> *<sup>K</sup>*5−<sup>4</sup> <sup>=</sup> <sup>−</sup>6(*EI*)<sup>2</sup> *L*2 2 , *<sup>K</sup>*4−<sup>6</sup> <sup>=</sup> *<sup>K</sup>*6−<sup>4</sup> <sup>=</sup> <sup>2</sup>(*EI*)<sup>2</sup> *<sup>L</sup>*<sup>2</sup> , *<sup>K</sup>*5−<sup>5</sup> <sup>=</sup> <sup>12</sup> (*EI*)<sup>2</sup> *L*3 2 + (*EI*)<sup>3</sup> *L*3 3 , *K*5−<sup>6</sup> = *K*6−<sup>5</sup> = 6 <sup>−</sup>(*EI*)<sup>2</sup> *L*2 2 + (*EI*)<sup>3</sup> *L*2 3 , *<sup>K</sup>*5−<sup>7</sup> <sup>=</sup> *<sup>K</sup>*7−<sup>5</sup> <sup>=</sup> <sup>−</sup>12(*EI*)<sup>3</sup> *L*3 3 , *<sup>K</sup>*5−<sup>8</sup> <sup>=</sup> *<sup>K</sup>*8−<sup>5</sup> <sup>=</sup> <sup>6</sup>(*EI*)<sup>3</sup> *L*2 3 , *K*6−<sup>6</sup> = 4 (*EI*)<sup>2</sup> *<sup>L</sup>*<sup>2</sup> <sup>+</sup> (*EI*)<sup>3</sup> *L*3 , *<sup>K</sup>*6−<sup>7</sup> <sup>=</sup> *<sup>K</sup>*7−<sup>6</sup> <sup>=</sup> <sup>−</sup>6(*EI*)<sup>3</sup> *L*2 3 , *<sup>K</sup>*6−<sup>8</sup> <sup>=</sup> *<sup>K</sup>*8−<sup>6</sup> <sup>=</sup> <sup>2</sup>(*EI*)<sup>3</sup> *<sup>L</sup>*<sup>3</sup> , *<sup>K</sup>*7−<sup>7</sup> <sup>=</sup> <sup>12</sup> (*EI*)<sup>3</sup> *L*3 3 + (*EI*)<sup>4</sup> *L*3 4 , *K*7−<sup>8</sup> = *K*8−<sup>7</sup> = 6 <sup>−</sup>(*EI*)<sup>3</sup> *L*2 3 + (*EI*)<sup>4</sup> *L*2 4 , *<sup>K</sup>*7−<sup>9</sup> <sup>=</sup> *<sup>K</sup>*9−<sup>7</sup> <sup>=</sup> <sup>−</sup>12(*EI*)<sup>4</sup> *L*3 4 , *<sup>K</sup>*7−<sup>10</sup> <sup>=</sup> *<sup>K</sup>*10−<sup>7</sup> <sup>=</sup> <sup>6</sup>(*EI*)<sup>4</sup> *L*2 4 , *K*8−<sup>8</sup> = 4 (*EI*)<sup>3</sup> *<sup>L</sup>*<sup>3</sup> <sup>+</sup> (*EI*)<sup>4</sup> *L*4 , *<sup>K</sup>*8−<sup>9</sup> <sup>=</sup> *<sup>K</sup>*9−<sup>8</sup> <sup>=</sup> <sup>−</sup>6(*EI*)<sup>4</sup> *L*2 4 , *<sup>K</sup>*8−<sup>10</sup> <sup>=</sup> *<sup>K</sup>*10−<sup>8</sup> <sup>=</sup> <sup>2</sup>(*EI*)<sup>4</sup> *<sup>L</sup>*<sup>4</sup> , *<sup>K</sup>*9−<sup>9</sup> <sup>=</sup> <sup>12</sup> (*EI*)<sup>4</sup> *L*3 4 + (*EI*)<sup>5</sup> *L*3 5 , *K*9−<sup>10</sup> = *K*10−<sup>9</sup> = 6 <sup>−</sup>(*EI*)<sup>4</sup> *L*2 4 + (*EI*)<sup>5</sup> *L*2 5 , *<sup>K</sup>*9−<sup>11</sup> <sup>=</sup> *<sup>K</sup>*11−<sup>9</sup> <sup>=</sup> <sup>−</sup>6(*L*5+2*L*6)(*EI*)<sup>5</sup> *L*3 <sup>5</sup>*L*<sup>6</sup> , *K*10−<sup>10</sup> = 4 (*EI*)<sup>4</sup> *<sup>L</sup>*<sup>4</sup> <sup>+</sup> (*EI*)<sup>5</sup> *L*5 , *<sup>K</sup>*10−<sup>11</sup> <sup>=</sup> *<sup>K</sup>*11−<sup>10</sup> <sup>=</sup> <sup>−</sup>2(*L*5+3*L*6)(*EI*)<sup>5</sup> *L*2 <sup>5</sup>*L*<sup>6</sup> , *<sup>K</sup>*11−<sup>11</sup> <sup>=</sup> *keq* <sup>+</sup> <sup>4</sup>(*L*<sup>2</sup> 5+3*L*5*L*6+3*L*<sup>2</sup> <sup>6</sup>)(*EI*)<sup>5</sup> *L*2 5*L*2 6

The electromechanical coupling matrix α is given by

$$\boldsymbol{\alpha} = \frac{bd\_{31}(t\_m + t\_p)}{s\_{11}} \begin{bmatrix} 0 & 1 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & -1 & 0 \end{bmatrix}^T$$

#### **References**


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

© 2020 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 (http://creativecommons.org/licenses/by/4.0/).
