Electromechanical finite element modelling for dynamic analysis of a cantilevered piezoelectric energy harvester with tip mass offset under base excitations

Recently, piezoelectric power harvesting devices have shown their viability for converting mechanical energy into useful electrical energy that can be used to power remote, smart wireless sensor systems. Their applications include monitoring of intelligent infrastructure, aircraft and environmental monitoring systems and biosensors for the human body. As a profound alternative energy source, a few published literature reviews [1–7] have discussed how the piezoelectric vibration-based power harvester applications from macro-to-micro levels benefit researchers who are investigating further novel techniques. Designing self-charging and lifespan energy harvesting devices requires technical knowledge and understanding of mechanical and electrical tuning techniques, material or physical properties, mechanical energy sources, power conditioning electronic circuits, sensor systems, and miniature fabrication. However, before conducting a manufacturing process for these smart devices in such conditions, initial mathematical models using analytical and numerical techniques become an important tool for analyzing the electromechanical piezoelectric vibration power harvesting response behaviour. Consequently, the development of new concepts has attracted attention from some researchers. Starting with the simplified analytical lumped parameter models of piezoelectric structures [8–10], the electrical equivalent systems representing the electromechanical piezoelectric component and standard circuit interface were developed. Considering the required accuracy and the necessity of modelling the cantilevered piezoelectric power harvester with tip mass under base excitation, the development of analytical approaches using Rayleigh–Ritz methods were explored in an effort to analyze various parametric case studies for vibration power harvesters [11–13]. Kim et al [14] investigated the vibrational piezoelectric bimorph power harvester with the effects of tip mass offset using the Rayleigh–Ritz method. Their interest and focus were to model the system responses of short- and open-circuit resonance frequency with variable load resistances using tip mass configuration. Moreover, Erturk [15] further developed comprehensive analytical approaches using assumed-mode methods for Euler–Bernoulli, Rayleigh and Timoshenko piezoelectric beams with axial deformations. It was proved that by increasing a certain number of modes in the choice of admissible trial functions, the approximation solution converges to the analytical solution, especially for the fundamental resonance frequency response. Erturk and Inman [16] have focused on researching the closed-form distributed parameter system of the Euler–Bernoulli bimorph beam for investigating frequency response functions (FRFs) under variable load resistance. Wickenheiser and Garcia [17] further developed the analytical modelling of piezoelectric bimorph beams with magnetic tip mass located near the fixed ferromagnetic structures in order to magnify the driving frequency below the fundamental frequency of the system. Moreover, Lumentut and Howard [18–20] explored the development of the piezoelectric bimorph beam with tip mass under both multidirectional excitation and the strain field effects of the transverse and initial longitudinal forms, in an effort to investigate a broad range of case studies using mathematical investigations of the orthonormalized electromechanical weak-form technique and closed-form boundary value method. In scrutinizing the system response of the single piezoelectric beam in light of the existing literature, the typical optimal power output mostly depends on the first mode. For this case, Lumentut et al [21] investigated the strategy of increasing power output and widening the multifrequency band using the orthonormalized multielectromechanical dynamic equations reduced from closed-form boundary value methods. Various parametric case studies have been developed using new electromechanical multimode FRFs of multiple piezoelectric bimorph beams connected in series, parallel and mixed series–parallel connections.

Some finite element literature concerns the use of piezoelectric components as smart material systems and has shown a variety of engineering applications, the most notable of which have been devoted to ultrasonic transducers and active control systems, either with feedback gain control or sensing and actuation response systems. However, for the current research, these piezoelectric finite element models give the relevancy and basis that allow the development of new numerical techniques for modelling the electromechanical power harvester, as presented in this paper. In the ultrasonic transducer application, Allik and Hughes [22] and Allik et al [23] first introduced mathematical models of the linear constitutive equations of the three-dimensional (3D) finite element piezoelectric vibration model using variational methods. Later, comprehensive mechanical discretized elements for analyzing an ultrasonic piezoelectric rod resonator [24] and the phenomenological electro-acoustic modelling for ultrasonic probes [25] were developed. Further application of the piezoelectric finite element has been developed for the active dynamic control system using feedback gain control. Tzou and Tseng [26] developed the active dynamic control matrix equations of the distributed piezoelectric components onto plate structures using the Hamiltonian principle and the Guyan reduction technique. Other published research works that analyze the controlled transient response [27–29] and transfer function [30] behaviour have been further investigated using the negative velocity feedback control system. For sensing and actuation response systems, the actuation system response, using applied voltage into the piezoelectric actuator, activated the piezoelectric sensor to create voltage output and deformation of the laminated composite structures [31–33]. More recently, Thomas et al [34] formulated the shunted electrical circuit connected to the piezoelectric patches using the finite element control system of the beam structure. Since the actual piezoelectric components are covered with electrode layers, the piezoelectric finite element literature ignores the electrode effect by assuming that the electrical voltages for each element and node are independent. Only a few of the published papers related to control finite element modelling have taken the electrode effect into consideration by assuming equal voltage distribution on each piezoelectric element [29, 31, 34]. However, the previous works are only based on mechanical discretization, without considering the electrical elements that are discretized due to the effect of electrode layers. Ignoring the electrical discretized elements from the electrode layer can affect the electromechanical FRF of the piezoelectric power harvester using the closed-circuit system. Also, there are no other previous research works that investigate the electromechanical finite element frequency responses using the electromechanical discretized elements under base excitation.

In earlier studies of the electromechanical finite element response of piezoelectric power harvesting, Lumentut et al [35] initially developed an energy harvesting plate structure with the segmented piezoelectric element based on the Love–Kirchhoff plate theory for investigating frequency and transient response behaviour. However, at an earlier stage, this model only considers the voltage output from the piezoelectric element without including the electrode effect with complete derivation of the electromechanical discretization related to the system response. DeMarqui et al [36] further presented the cantilevered plate with distributed piezoelectric layers by including the electrode effect, as validated by Erturk and Inman [16] in the existing literature of the piezoelectric bimorph with tip mass. However, their proposed finite element model does not take into account the important physical issues of the electromechanical orthonormalized global matrix and scalar forms, and the mathematical techniques of electromechanical discretization. Aladwani et al [37] further investigated the dynamic magnifier using a spring system attached to the base structure of the cantilevered piezoelectric bimorph beam with tip mass. The system was used to model the optimized power harvesting frequency response behaviour. However, their proposed work neglects the electrode effect of the piezoelectric component for formulating the electromechanical discretized elements, and also ignores the offset effect of the tip mass measured between the tip mass centroid and the end of the bimorph. Further research investigating piezoelectric power harvesters has been conducted using solid finite element-based ANSYS multiphysics software. Zhu et al [38] investigated the prediction of power output of FRFs under the resistive shunt circuit. Moreover, Yang and Tang [39] developed the electrical equivalent circuit model using SPICE software, whose circuit parameters were adopted from the ANSYS simulation for predicting the power output. However, although ANSYS software provides capabilities of 3D, coupled field solid, and plane elements, it is not able to model the laminated beam elements. For the standard level of the energy harvesters, mesh generation using laminated beam elements is sufficient to model the laminated piezoelectric beam structures (e.g., unimorph, bimorph and multimorph).

In this paper, new numerical techniques of the electromechanical finite element vibration power harvester using the Euler–Bernoulli unimorph beam model are presented and validated with experimental studies. The numerical models proposed here provide new technical contributions for analyzing the interrelationship between coupled field equations, electromechanical discretization, kinematic motions, Lagrangian electromechanical dynamic equations and electromechanical finite element equations. These technical contributions provide clear insight into the new electromechanical finite element equations for power harvesting applications, which can be highlighted in the following points:

• In considering the effect of electrode layers, we introduce the new electromechanical discretization, consisting of mechanical and electrical discretized elements, to model the electromechanical finite element piezoelectric structure. The electromechanical discretized elements provide an important basis for formulating the FRFs using new electromechanical finite element dynamic equations.
• The elemental system dynamics of the piezoelectric unimorph beam, including tip mass offset, are reduced from the kinematic equations. It is found that the effect of the offset distance, measured from the tip mass centroid and the end of the unimorph beam, can result in extra terms in the mass matrix and the generalized input dynamic force vector. The benefit of applying the tip mass offset, located after the end of the piezoelectric length, is that it can provide a very practical MEMS power harvester design for any level of scalability. It can also avoid the use of material around the end of the piezoelectric length, which can readily be damaged because of its brittle nature.
• The established, new, electromechanical finite element approach can be used conveniently to model different levels of scalability, including the segmented piezoelectric element coupled with substructure and micro-electro-mechanical system (MEMS) devices.
• The proposed numerical techniques also provide a complement to study the laminated piezoelectric beam with electrode layers in different applications.
• Lagrangian electromechanical dynamic equations are developed to formulate the constitutive elemental matrices of the electromechanical finite element equations. The orthonormalized global matrices and global scalar forms can be further reduced to formulate new multimode electromechanical FRFs of the displacement, velocity, voltage, current and power harvesting, including the optimal power harvester. The normalized strain node and eigenmode shapes are formulated using discretized elements for convergence studies. Finally, the parametric case studies with resistive shunt circuits, based on the proposed novel electromechanical finite element and experimental validations, are also presented and discussed.

Constitutive linear piezoelectric beam equations can be formulated using tensor electrical enthalpy concepts, which are based on the continuum thermodynamics that can be condensed using Voigt's notation and Einstein's summation convention [19, 40–42]. The electrical enthalpy of typical smart materials with lead zirconate titanate (PZT) can be formulated under adiabatic and isothermal processes as

$$\begin{eqnarray}&&\hat{H}\left( {{{\boldsymbol{S}} }_{1}},{{{\boldsymbol{E}} }_{3}} \right)=\frac{1}{2}\bar{c}_{11}^{E}{{{\boldsymbol{S}} }_{1}}{{{\boldsymbol{S}} }_{1}}-{{e}_{31}}{{{\boldsymbol{E}} }_{3}}{{{\boldsymbol{S}} }_{1}}-\frac{1}{2}\varepsilon _{33}^{S}{{{\boldsymbol{E}} }_{3}}{{{\boldsymbol{E}} }_{3}}.\end{eqnarray} \tag{ 1 }$$

Notations of the piezoelectric parameters used here are stated in accordance with IEEE standards [43]. Here, the parameters $\bar{c}_{11}^{E}$, ${{e}_{{\rm 31}}}$, $\;\varepsilon _{33}^{S}$, ${{{\boldsymbol{E}} }_{3}}$, ${{{\boldsymbol{T}} }_{1\,}}$ and ${{{\boldsymbol{S}} }_{1}}$ represent the piezoelectric elastic stiffness at constant electric field, piezoelectric coefficient, permittivity under constant strain, electric field, stress and strain, respectively. The parameter $\varepsilon _{33}^{S}=\varepsilon _{33}^{T}-{{e}_{31}}{{d}_{31}}$ or $\varepsilon _{33}^{S}=\varepsilon _{33}^{T}-d_{31}^{2}\;\bar{c}_{11}^{E}$ is given where $\bar{c}_{11}^{E}\;=\;1/s_{11}^{E}$, $\varepsilon _{33}^{T}$ is the permittivity at constant stress and $s_{11}^{E}$ is the elastic compliance at constant electric field. Moreover, the common piezoelectric constant produced from the manufacturing company is in the form ${{d}_{31}}$, where this can be modified by multiplying the plane stress-based elastic stiffness at constant electric field to give ${{e}_{31}}\;=\;{{d}_{31}}\bar{c}_{11}^{E}$. The electric displacement of the piezoelectric element, ${{D}_{3}}$, can be obtained by differentiating equation (1) with respect to E ₃. This electric displacement of piezoelectricity can be formulated as

$$\begin{eqnarray}&&{{{\boldsymbol{D}} }_{3}}=-\frac{\partial \hat{H}}{\partial {{{\boldsymbol{E}} }_{3}}}={{e}_{31}}{{{\boldsymbol{S}} }_{1}}+\varepsilon _{33}^{S}{{{\boldsymbol{E}} }_{3}}.\end{eqnarray} \tag{ 2 }$$

Note that for the purpose of convenience in the forthcoming notations, the superscripts 1 and 2 that are used in the parameters of stress, strain, elastic stiffness, density and cross sectional area refer to the substructure and piezoelectric layers, respectively. Since the piezoelectric unimorph consists of active and inactive layers, the plane stress field from the inactive or substructure layer can also be formulated as

$$\begin{eqnarray}&&{\boldsymbol{T}} _{1}^{\left( 1 \right)}=\bar{c}_{11}^{\left( 1 \right)}{\boldsymbol{S}} _{1}^{\left( 1 \right)}.\end{eqnarray} \tag{ 3 }$$

The piezoelectric stress field can also be formulated by differentiating equation (1) with respect to S _1, to give

$$\begin{eqnarray}&&{\boldsymbol{T}} _{1}^{\left( 2 \right)}=\frac{\partial \hat{H}}{\partial {\boldsymbol{S}} _{1}^{\left( 2 \right)}}=\bar{c}_{11}^{\left( 2 \right)}{\boldsymbol{S}} _{1}^{\left( 2 \right)}-{{e}_{31}}{{{\boldsymbol{E}} }_{3}}.\end{eqnarray} \tag{ 4 }$$

Note that alternative derivations of equations (2) and (4) can also be obtained by using the enthalpy equation of state in terms of continuum thermopiezoelectricity, Maxwell's relations and Legendre's transformation [41, 44].

In this section, the kinematic motions for a continuous geometrical beam, as shown in figure 1, are derived from the different configurations at a series of points, resulting in the deformations and velocity of the body. However, the kinematics for the tip mass offset, whose motions are solely affected by the tip of the beam's transverse rectilinear $[\dot{w}_{abs}^{tip}\left( L,t \right)]$ and angular $[\dot{\theta }\left( L,t \right)]$ velocities, is a rigid-body motion with different configurations that result in extra terms in the kinematic equations. In reviewing previous studies, the most recently published papers on MEMS power harvesters [45–50] included the tip mass offset without scrutinizing its effect on the kinematic equations. The effect of the offset distance between the tip mass centroid and the end of the beam can effectively contribute to the mass matrix and the generalized input dynamic forces, as shown in section 4.3. It is important to note that the aim of formulating the kinematic equations is to formulate the Lagrangian functional energy forms that consist of the kinetic, potential and electrical energies of the piezoelectric unimorph, with the tip mass as given in the forthcoming section.

Zoom In Zoom Out Reset image size

In figure 1, the undeformed beam under base vector ${{w}_{base}}\left( t \right)$ moves from the fixed reference frame of oXZ to the initial reference frame of $o^{\prime} X^{\prime} Z$. Since the base vector, as the input excitation, moves from point $o$ to $o^{\prime}$ at the designated reference frame, the position of point $p$ also moves to point $p^{\prime}$ with the same magnitude. As a result, the base vector creates the relative transverse deformation, $w\left( x,t \right)$, measured from point $p^{\prime}$ moving to the final point, $p^{\prime\prime}$. The absolute displacement, ${{w}_{abs}}\left( x,t \right)$, can also be obtained from the fixed reference frame of oXZ to the final position. Note that the difference between the absolute displacement and the base vector defines the relative deformation. Moreover, the beam also carries an arbitrary tip mass offset whose kinematic motions must be considered. The position vectors of the tip mass start at point d to the final differential element at point $m^{\prime\prime}$, where the resultant vectors at points $o-d^{\prime\prime} -g^{\prime\prime} -m^{\prime\prime}$ can be defined. The point of attachment, $d^{\prime\prime}$ (deformation point at the end of the beam), to the tip mass centroid at point $g^{\prime\prime}$ can be viewed as the offset of the tip mass. The origin of tip mass centroid can be determined from the differential element, using the extra position vector from points $g^{\prime\prime}$ to $m^{\prime\prime}$. As mentioned previously, since the tip mass was assumed to be a rigid body, its relative motions depend on the tip of the beam's motion at point $d^{\prime\prime}$.

The position vector, ${{{\boldsymbol{R}} }^{pp^{\prime\prime} }}$, for the elemental beam can be defined as

$$\begin{eqnarray}&&{{{\boldsymbol{R}} }^{pp^{\prime\prime} }}\left( x,z,t \right)={{{\boldsymbol{R}} }^{op^{\prime\prime} }}-{{{\boldsymbol{R}} }^{op}},\end{eqnarray} \tag{ 5 }$$

where the position vectors, ${{{\boldsymbol{R}} }^{op^{\prime\prime} }}$ and ${{{\boldsymbol{R}} }^{op}}$, can be written as

$$\begin{eqnarray}\begin{array}{lllllllllllllll} {{{\boldsymbol{R}} }^{op^{\prime\prime} }} & = & \left( x-z\;{\rm sin} \;\theta \left( x,t \right) \right){{{\boldsymbol{e}} }_{1}} \\ {} & {} & +\left( {{w}_{abs}}\left( x,t \right)+z\;{\rm cos} \;\theta \left( x,t \right) \right){{{\boldsymbol{e}} }_{3}}, \\ {{{\boldsymbol{R}} }^{op}} & = & x{{{\boldsymbol{e}} }_{1}}+z{{{\boldsymbol{e}} }_{3}}. \\ \end{array}\end{eqnarray} \tag{ 6 }$$

In simplification of the linear assumption, the small angle, $\theta \left( x,t \right)=\partial w\left( x,t \right)/\partial x$, can be obtained after applying Taylor's theorem to give ${\rm sin} \;\theta \left( x,t \right)\approx \theta \left( x,t \right)$ and ${\rm cos} \;\theta \left( x,t \right)\approx 1$. As the initial position vector of the elemental beam, ${{{\boldsymbol{R}} }^{op}}$, is defined, the position vector, ${{{\boldsymbol{R}} }^{pp^{\prime\prime} }}$, from equation (5) can be differentiated with respect to time, which gives

$$\begin{eqnarray}\begin{array}{lllllllllllllll} {{{{\boldsymbol{\dot{R}}} }}^{pp^{\prime\prime} }}\left( x,z,t \right) & = & {{{{\boldsymbol{\dot{R}}} }}^{op^{\prime\prime} }}-{{{{\boldsymbol{\dot{R}}} }}^{op}} \\ {} & = & -\dot{\theta }\left( x,t \right){{{\boldsymbol{e}} }_{2}}\times z{{{\boldsymbol{e}} }_{3}} \\ {} & {} & +\left( {{{\dot{w}}}_{base}}\left( t \right)+\dot{w}\left( x,t \right) \right){{{\boldsymbol{e}} }_{3}}. \\ \end{array}\end{eqnarray} \tag{ 7 }$$

Note that the absolute velocity vectors for the elemental beam and tip mass can be written as

$$\begin{eqnarray}\begin{array}{lllllllllllllll} {{{{\boldsymbol{\dot{w}}} }}_{abs}}\left( x,t \right) & = & \left( {{{\dot{w}}}_{base}}\left( t \right)+\dot{w}\left( x,t \right) \right){{{\boldsymbol{e}} }_{3}} \\ {\boldsymbol{\dot{w}}} _{abs}^{tip}\left( L,t \right) & = & \left( {{{\dot{w}}}_{base}}\left( t \right)+\dot{w}\left( L,t \right) \right){{{\boldsymbol{e}} }_{3}}. \\ \end{array}\end{eqnarray} \tag{ 8 }$$

For the elemental tip mass, the position vector, ${{{\boldsymbol{R}} }^{m{{m}^{\prime\prime }}}}$, can be formulated as

$$\begin{eqnarray}\begin{array}{lllllllllllllll} {{{\boldsymbol{R}} }^{m{{m}^{\prime\prime }}}}\left( L,z,t \right) & = & {{{\boldsymbol{R}} }^{o{{d}^{\prime\prime }}}}+{{{\boldsymbol{R}} }^{{{d}^{\prime\prime }}{{g}^{\prime\prime }}}}+{{{\boldsymbol{R}} }^{{{g}^{\prime\prime }}{{m}^{\prime\prime }}}}-{{{\boldsymbol{R}} }^{od}}-{{{\boldsymbol{R}} }^{dg}}-{{{\boldsymbol{R}} }^{gm}}, \\ \end{array}\end{eqnarray} \tag{ 9 }$$

where each position vector can be formulated as

$$\begin{eqnarray}\begin{array}{lllllllllllllll} {{{\boldsymbol{R}} }^{od}} & = & L{{{\boldsymbol{e}} }_{1}}, \\ {{{\boldsymbol{R}} }^{dg}} & = & {{x}_{c}}{{{\boldsymbol{e}} }_{1}}+{{z}_{c}}{{{\boldsymbol{e}} }_{3}}, \\ {{{\boldsymbol{R}} }^{gm}} & = & {{x}_{m}}{{{\boldsymbol{e}} }_{1}}+{{z}_{m}}{{{\boldsymbol{e}} }_{3}}, \\ {{{\boldsymbol{R}} }^{o{{d}^{\prime\prime }}}} & = & L{{{\boldsymbol{e}} }_{1}}+\left( {{w}_{base}}\left( t \right)+w\left( L,t \right) \right){{{\boldsymbol{e}} }_{3}} \\ {{{\boldsymbol{R}} }^{{{d}^{\prime\prime }}{{g}^{\prime\prime }}}} & = & \left( {{x}_{c}}{\rm cos} \theta \left( L,t \right)-{{z}_{c}}{\rm sin} \theta \left( L,t \right) \right){{{\boldsymbol{e}} }_{1}} \\ {} & {} & +\left( {{z}_{c}}{\rm cos} \theta \left( L,t \right)+{{x}_{c}}{\rm sin} \theta \left( L,t \right) \right){{{\boldsymbol{e}} }_{3}} \\ {{{\boldsymbol{R}} }^{{{g}^{\prime\prime }}{{m}^{\prime\prime }}}} & = & \left( {{x}_{m}}\;{\rm cos} \;\theta \left( L,t \right)-{{z}_{m}}\;{\rm sin} \;\theta \left( L,t \right) \right){{{\boldsymbol{e}} }_{1}} \\ {} & {} & +\left( {{z}_{m}}\;{\rm cos} \;\theta \left( L,t \right)+{{x}_{m}}\;{\rm sin} \;\theta \left( L,t \right) \right){{{\boldsymbol{e}} }_{3}}. \\ \end{array}\end{eqnarray} \tag{ 10 }$$

Again, the linear system for the small angle, $\theta \left( L,t \right)=\partial w\left( L,t \right)/\partial x$, can be made using Taylor's theorem to give ${\rm sin} \;\theta \left( L,t \right)\approx \theta \left( L,t \right)$ and ${\rm cos} \;\theta \left( L,t \right)\approx 1$. In this case, the velocity of the elemental tip mass can be obtained by differentiating equation (9) with respect to time, to yield

$$\begin{eqnarray}\begin{array}{lllllllllllllll} {{{{\boldsymbol{\dot{R}}} }}^{m{{m}^{\prime\prime }}}}\left( L,z,t \right) & = & {{{{\boldsymbol{\dot{R}}} }}^{o{{d}^{\prime\prime }}}}+{{{{\boldsymbol{\dot{R}}} }}^{{{d}^{\prime\prime }}{{g}^{\prime\prime }}}}+{{{{\boldsymbol{\dot{R}}} }}^{{{g}^{\prime\prime }}{{m}^{\prime\prime }}}}-{{{{\boldsymbol{\dot{R}}} }}^{od}}-{{{{\boldsymbol{\dot{R}}} }}^{dg}}-{{{{\boldsymbol{\dot{R}}} }}^{gm}} \\ {} & = & \left( {{{\dot{w}}}_{base}}\left( t \right)+\dot{w}\left( L,t \right) \right){{{\boldsymbol{e}} }_{3}} \\ {} & {} & -\dot{\theta }\left( L,t \right){{{\boldsymbol{e}} }_{2}}\left( {{z}_{c}}{{{\boldsymbol{e}} }_{3}}+{{x}_{c}}{{{\boldsymbol{e}} }_{1}} \right) \\ {} & {} & -\dot{\theta }\left( L,t \right){{{\boldsymbol{e}} }_{2}}\left( {{z}_{m}}{{{\boldsymbol{e}} }_{3}}+{{x}_{m}}{{{\boldsymbol{e}} }_{1}} \right). \\ \end{array}\end{eqnarray} \tag{ 11 }$$

The position vector, ${{{\boldsymbol{R}} }^{{{p}^{\prime }}{{p}^{\prime\prime }}}}$, can be specified as the relative displacement with respect to the moving base support from reference frame $oXZ$ to $o^{\prime} X^{\prime} Z^{\prime}$, as

$$\begin{eqnarray}\begin{array}{lllllllllllllll} {{{\boldsymbol{R}} }^{p^{\prime} {{p}^{\prime\prime }}}}\left( x,z,t \right) & = & {{{\boldsymbol{R}} }^{p{{p}^{\prime\prime }}}}-{{{\boldsymbol{R}} }^{pp^{\prime} }} \\ {} & = & -z\theta \left( x,t \right){{{\boldsymbol{e}} }_{1}}+w\left( x,t \right){{{\boldsymbol{e}} }_{3}}. \\ \end{array}\end{eqnarray} \tag{ 12 }$$

The elemental beam strain can be obtained by differentiating ${{{\boldsymbol{R}} }^{{{p}^{\prime }}{{p}^{\prime\prime }}}}$ with respect to x, giving the typical Euler–Bernoulli beam theory as

$$\begin{eqnarray}&&{{S}_{1}}=\frac{\partial ({{{\boldsymbol{R}} }^{p^{\prime} {{p}^{\prime\prime }}}}\cdot {{{\boldsymbol{e}} }_{1}})}{\partial x}=-z\frac{{{\partial }^{2}}w\left( x,t \right)}{\partial {{x}^{2}}}.\end{eqnarray} \tag{ 13 }$$

Since the piezoelectric unimorph beam is viewed as the asymmetric structure with different material properties and geometrical dimensions, the asymmetric neutral axis [51] must be determined as presented in appendix A.

In common practice, the surfaces of actual or physical piezoelectric material currently on the market are covered evenly with thin conducting electrodes, resulting in single voltage output during dynamic motions. As shown in figure 2(a), we chose the physical structure of the unimorph beam with an arbitrary tip mass offset, which consists of a piezoelectric material with a conducting electrode layer and a substructure. However, in figure 2(b), the electromechanical discretization, consisting of mechanical and electrical discretized finite elements, is completely different from the physical structure. In this situation, the electromechanical finite element discretization of the piezoelectric structure should be analyzed using the generalized multi-output electrical current from each element connected with an electrical parallel connection. Since the polarization is proportional to the stress field and the stress field is also proportional to the strain field [40, 44, 52], the polarization, P , and electric field, E , directions are chosen in the three-direction (thickness) along the z-axis.

Zoom In Zoom Out Reset image size

4.1. Geometrical analysis of an electromechanical element

The first-order Hermite interpolation of the cubic displacement function is chosen here to formulate the unknown four degrees of freedom for the two-node unimorph beam element [53–55]. The cubic displacement function must satisfy the beam theory continuity conditions for the translation and rotational displacements for the general location of the elements. It is important to note here that the kinematic equations given in section 3 can be used to formulate the local element near the tip mass, as shown in figure 3, since the translation and angular velocities of the tip massʼs rigid body depend on the end of the unimorph beam of length, giving particular insight to formulating the complete electromechanical finite element equations.

Zoom In Zoom Out Reset image size

The first-order Hermite interpolation in terms of the polynomial, ${{\Phi }^{\left( 1 \right)}}\left( x \right)$, and the unknown nodal displacement, ${{u}_{i}}\left( t \right)$, for the elemental beam can be formulated as

$$\begin{eqnarray}&&w\left( x,t \right)=\sum \limits_{i=1}^{2} \sum \limits_{k=0}^{1} \Phi _{ki}^{\left( 1 \right)}\left( x \right)u_{i}^{\left( k \right)}\left( t \right).\end{eqnarray} \tag{ 14 }$$

The property of Hermite polynomials for modelling the translation and rotation at the nodes can be formulated in terms of the l-derivatives to give

$$\begin{eqnarray}\begin{array}{lllllllllllllll} \frac{{{{\rm d}}^{l}}\Phi _{ki}^{\left( 1 \right)}}{{\rm d}{{x}^{l}}}\left( {{x}_{n}} \right) & = & {{\delta }_{in}}{{\delta }_{kl}}{\rm for}\;i,n=1,{\mkern 1mu} 2\;{\rm and} \\ k,l & = & 0,1 \\ \end{array},\end{eqnarray} \tag{ 15 }$$

where x_n is the location at the nth node and ${{\delta }_{pq}}$ is the Kronecker delta with the property

$$\begin{eqnarray}&&{{\delta }_{pq}}=\left\{ \begin{array}{ccccccccccccccc} 0\quad if\;p\ne q \\ 1\quad if\;p=q \\ \end{array} \right..\end{eqnarray} \tag{ 16 }$$

Corresponding to equations (15) and (16), the lth derivatives of w(x, t) can be written as

$$\begin{eqnarray}&&\frac{{{{\rm d}}^{l}}{{w}_{n}}}{{\rm d}{{x}^{l}}}\left( x,t \right)=\sum \limits_{i=1}^{2} \sum \limits_{k=0}^{1} \frac{{{{\rm d}}^{l}}\Phi _{ki}^{\left( 1 \right)}}{{\rm d}{{x}^{l}}}\left( {{x}_{n}} \right)u_{i}^{\left( k \right)}\left( t \right),\end{eqnarray} \tag{ 17 }$$

or they can be formulated in terms of the unknown nodal displacement forms as

$$\begin{eqnarray}&&\frac{{{{\rm d}}^{l}}{{w}_{n}}}{{\rm d}{{x}^{l}}}\left( t \right)=\sum \limits_{i=1}^{2} \sum \limits_{k=0}^{1} {{\delta }_{in}}{{\delta }_{kl}}u_{i}^{\left( k \right)}\left( t \right)=u_{n}^{\left( l \right)}\left( t \right).\end{eqnarray} \tag{ 18 }$$

Note that equation (18) clearly shows certain physical meaning after applying equations (15) and (16), giving the nodal displacement forms as a function of time. Corresponding to index notation from equation (15), equation (14) can be modified by applying equation (18) to give

$$\begin{eqnarray}&&w\left( x,t \right)=\sum \limits_{i=1}^{2} \sum \limits_{k=0}^{1} \Phi _{ki}^{\left( 1 \right)}\left( x \right)\frac{{{{\rm d}}^{k}}{{w}_{i}}}{{\rm d}{{x}^{k}}}\left( t \right).\end{eqnarray} \tag{ 19 }$$

To find four terms of first-order Hermite polynomials, the property given from equation (15) can be used to identify four elemental boundary conditions for each degree of freedom at the nodes. In such conditions, the cubic equations with four unknown coefficients for each degree of freedom at the nodes can be used, and the results of the Hermite polynomials can be formulated after simplifying to give

$$\begin{eqnarray}\begin{array}{lllllllllllllll} \Phi _{01}^{\left( 1 \right)}\left( x \right) & = & 1-3{{\left( \frac{x-{{x}_{e}}}{{{L}_{e}}} \right)}^{2}}+2{{\left( \frac{x-{{x}_{e}}}{{{L}_{e}}} \right)}^{3}}, \\ \Phi _{11}^{\left( 1 \right)}\left( x \right) & = & \left( x-{{x}_{e}} \right){{\left( 1-\frac{x-{{x}_{e}}}{{{L}_{e}}} \right)}^{2}} \\ \Phi _{02}^{\left( 1 \right)}\left( x \right) & = & 3{{\left( \frac{x-{{x}_{e}}}{{{L}_{e}}} \right)}^{2}}-2{{\left( \frac{x-{{x}_{e}}}{{{L}_{e}}} \right)}^{3}}, \\ \Phi _{12}^{\left( 1 \right)}\left( x \right) & = & \left( x-{{x}_{e}} \right)\left[ {{\left( \frac{x-{{x}_{e}}}{{{L}_{e}}} \right)}^{2}}-\frac{x-{{x}_{e}}}{{{L}_{e}}} \right] \\ {{L}_{e}} & = & {{x}_{e+1}}-{{x}_{e}}. \\ \end{array}\end{eqnarray} \tag{ 20 }$$

Corresponding to equation (14), the simplified first-order Hermite interpolation function of the unimorph beam can be modified into a matrix form to give

$$\begin{eqnarray}&&{\boldsymbol{w}} \left( x,t \right)={{{\boldsymbol{ \Phi }} }^{e}}\left( x \right){{{\boldsymbol{u}} }^{e}}\left( t \right),\end{eqnarray} \tag{ 21 }$$

where

$$\begin{eqnarray}\begin{array}{lllllllllllllll} {{{\boldsymbol{ \Phi }} }^{e}}\left( x \right) & = & \left[ {{\Phi }_{1}}\left( x \right){{\Phi }_{2}}\left( x \right){{\Phi }_{3}}\left( x \right){{\Phi }_{4}}\left( x \right) \right] \\ {{{\boldsymbol{u}} }^{e}}\left( t \right) & = & {{\left[ {{u}_{1}}{{u}_{2}}{{u}_{3}}{{u}_{4}} \right]}^{T}}. \\ \end{array}\end{eqnarray} \tag{ 22 }$$

Each parameter from equation (22) can be written in terms of the elemental displacement vector, u, and shape function, $\;\Phi$, for each node as

$$\begin{eqnarray}\begin{array}{lllllllllllllll} {{u}_{1}}\left( t \right) & = & u_{1}^{\left( 0 \right)}\left( t \right)={{w}_{1}}\left( t \right), \\ {{u}_{2}}\left( t \right) & = & u_{1}^{\left( 1 \right)}\left( t \right)={{\theta }_{1}}\left( t \right), \\ {{u}_{3}}\left( t \right) & = & u_{2}^{\left( 0 \right)}\left( t \right)={{w}_{2}}\left( t \right), \\ {{u}_{4}}\left( t \right) & = & u_{2}^{\left( 1 \right)}\left( t \right)={{\theta }_{2}}\left( t \right) \\ {{\Phi }_{1}}\left( x \right) & = & \Phi _{01}^{\left( 1 \right)}\left( x \right), \\ {{\Phi }_{2}}\left( x \right) & = & \Phi _{11}^{\left( 1 \right)}\left( x \right), \\ {{\Phi }_{3}}\left( x \right) & = & \Phi _{02}^{\left( 1 \right)}\left( x \right), \\ {{\Phi }_{4}}\left( x \right) & = & \Phi _{12}^{\left( 1 \right)}\left( x \right). \\ \end{array}\end{eqnarray} \tag{ 23 }$$

Corresponding to equation (13), the strain–displacement relationship, in terms of the differential form of the shape function of the piezoelectric unimorph beam, can be formulated by substituting equation (21) to give

$$\begin{eqnarray}&&{\boldsymbol{S}} \left( x,t \right)=-z{{{\boldsymbol{ \Psi }} }^{e}}\left( x \right){{{\boldsymbol{u}} }^{e}}\left( t \right),\end{eqnarray} \tag{ 24 }$$

where the differential form of the shape function of the strain–displacement relationship can be formulated as

$$\begin{eqnarray}\begin{array}{lllllllllllllll} {{{\boldsymbol{ \Psi }} }^{e}}\left( x \right) & = & \frac{{{{\rm d}}^{2}}{{{\boldsymbol{ \Phi }} }^{e}}\left( x \right)}{{\rm d}{{x}^{2}}} \\ {} & = & \left[ \begin{array}{ccccccccccccccc} {{\Psi }_{1}}\left( x \right) & {{\Psi }_{2}}\left( x \right) & {{\Psi }_{3}}\left( x \right) & {{\Psi }_{4}}\left( x \right) \\ \end{array} \right], \\ \end{array}\end{eqnarray} \tag{ 25 }$$

where

$$\begin{eqnarray}\begin{array}{lllllllllllllll} {{\Psi }_{1}}\left( x \right) & = & -\frac{6}{L_{e}^{2}}+\frac{12\left( x-{{x}_{e}} \right)}{L_{e}^{3}}, \\ {{\Psi }_{2}}\left( x \right) & = & \frac{6\left( x-{{x}_{e}} \right)}{L_{e}^{2}}-\frac{4}{{{L}_{e}}} \\ {{\Psi }_{3}}\left( x \right) & = & -\frac{12\left( x-{{x}_{e}} \right)}{L_{e}^{3}}+\frac{6}{L_{e}^{2}}, \\ {{\Psi }_{4}}\left( x \right) & = & \frac{6\left( x-{{x}_{e}} \right)}{L_{e}^{2}}-\frac{2}{{{L}_{e}}}. \\ \end{array}\end{eqnarray} \tag{ 26 }$$

The elemental discretized electric field, E , induced by the strain field due to input ambient vibration, can create a polarization in the piezoelectric material in the z-direction along its thickness, generating the electrical voltages. The electrical field is a function of the electrical potential, with a negative gradient operator given by

$$\begin{eqnarray}&&{{{\boldsymbol{E}} }_{3}}=-\nabla {{\phi }^{e}}\left( z,t \right),\end{eqnarray} \tag{ 27 }$$

where ${{\phi }^{e}}\left( z,t \right)={{\vartheta }^{e}}\left( z \right){{{\boldsymbol{v}} }^{e}}\left( t \right)$ is the electrical potential with linear assumption, the parameter ${{\vartheta }^{e}}\left( z \right)=(z-{{z}_{n}}+{{h}_{p}})/{{h}_{p}}$ is the shape function over the interval ${{z}_{n}}-{{h}_{p}}\leqslant z\leqslant {{z}_{n}}$, and z_n indicates the distance from the asymmetric neutral axis, as given in appendix A. The symbol $\nabla$ is a gradient operator for the first derivative of the shape function with respect to the thickness direction, giving ${{{\boldsymbol{ \Omega }} }^{\left( e \right)}}\left( z \right)\;=\;d{{\vartheta }^{e}}\left( z \right)/dz=\;\;1/{{h}_{p}}$. In this case, the discretized electric field can be reformulated as

$$\begin{eqnarray}&&{{{\boldsymbol{E}} }_{3}}=-{{{\boldsymbol{ \Omega }} }^{e}}\left( z \right){{{\boldsymbol{v}} }^{\left( e \right)}}\left( t \right).\end{eqnarray} \tag{ 28 }$$

The expression of the stress field in partial differential shape function form can be formulated for both piezoelectric and substructure elements by substituting equations (24) and (28) into equations (3) and (4), to give

$$\begin{eqnarray}\begin{array}{lllllllllllllll} {\boldsymbol{T}} _{1}^{\left( 2 \right)} & = & -z\bar{c}_{11}^{\left( 2 \right)}{{{\boldsymbol{ \Psi }} }^{e}}\left( x \right){{{\boldsymbol{u}} }^{e}}\left( t \right)+{{e}_{31}}{{{\boldsymbol{ \Omega }} }^{e}}\left( z \right){{{\boldsymbol{v}} }^{e}}\left( t \right) \\ {\boldsymbol{T}} _{1}^{\left( 1 \right)} & = & -z\bar{c}_{11}^{\left( 1 \right)}{{{\boldsymbol{ \Psi }} }^{e}}\left( x \right){{{\boldsymbol{u}} }^{e}}\left( t \right). \\ \end{array}\end{eqnarray} \tag{ 29 }$$

The electric displacement vector of the piezoelectric component can be formulated by substituting equations (24) and (28) into equation (2) to give

$$\begin{eqnarray}&&{{{\boldsymbol{D}} }_{3}}=z{{e}_{31}}{{{\boldsymbol{ \Psi }} }^{e}}\left( x \right){{{\boldsymbol{u}} }^{e}}\left( t \right)-\varepsilon _{33}^{S}{{{\boldsymbol{ \Omega }} }^{e}}\left( z \right){{{\boldsymbol{v}} }^{e}}\left( t \right).\end{eqnarray} \tag{ 30 }$$

4.2. Lagrangian electromechanical finite element equation

Constitutive electromechanical finite element dynamic equations are developed for modelling the elemental piezoelectric unimorph with the tip mass offset under input base transverse excitation. The functional energy forms using the extended Lagrangian formulation consist of kinetic (KE), potential (PE) and electrical (WE) energies, and nonconservative work (WF). The extended Lagrangian equation can be written as

$$\begin{eqnarray}&&\mathop \sum \limits_{s=1} \left[ \frac{{\rm d}}{{\rm d}t}\frac{\partial \Pi }{\partial {{{\boldsymbol{ \eta }} }_{s}}}-\frac{\partial \Pi }{\partial {{{\boldsymbol{ \eta }} }_{s}}}-{{{\boldsymbol{F}} }_{s}} \right]=0,\end{eqnarray} \tag{ 31 }$$

where $\Pi =KE\;-PE\;+\;WE$, ${{{\boldsymbol{ \eta }} }_{s}}\in \{{\boldsymbol{\dot{u}}} \left( t \right),{\boldsymbol{u}} \left( t \right),{\boldsymbol{v}} \left( t \right)\}$ and ${{F}_{s}}\in \{{{{\boldsymbol{\ddot{w}}} }_{base}}\left( t \right),{\boldsymbol{q}} \left( t \right)\}$ or $\{\sum s=1 {{F}_{s}}\}=\{\frac{\partial WF}{\partial {{{\boldsymbol{u}} }^{e}}\left( t \right)}+\frac{\partial WF}{\partial {{{\boldsymbol{v}} }^{e}}\left( t \right)}\}$.

The kinetic energy can be formulated based on the mass density of the piezoelectric unimorph and the density of tip mass, as

$$\begin{eqnarray}\begin{array}{lllllllllllllll} KE & = & \frac{1}{2}\int \limits_{{{x}_{e}}}^{{{x}_{e+1}}} \mathop \int \limits_{{{A}^{\left( 1 \right)}}} {{z}^{2}}{{\rho }^{\left( 1 \right)}}{\boldsymbol{\dot{ \theta }}} {{\left( x,t \right)}^{T}}{\boldsymbol{\dot{ \theta }}} \left( x,t \right){\rm d}{{A}^{\left( 1 \right)}}{\rm d}x \\ {} & {} & +\frac{1}{2}\int \limits_{{{x}_{e}}}^{{{x}_{e+1}}} \mathop \int \limits_{{{A}^{\left( 2 \right)}}} {{z}^{2}}{{\rho }^{\left( 2 \right)}}{\boldsymbol{\dot{ \theta }}} {{\left( x,t \right)}^{T}}{\boldsymbol{\dot{ \theta }}} \left( x,t \right){\rm d}{{A}^{\left( 2 \right)}}{\rm d}x \\ {} & {} & +\frac{1}{2}\int \limits_{{{x}_{e}}}^{{{x}_{e+1}}} \mathop \int \limits_{{{A}^{\left( 1 \right)}}} {{\rho }^{\left( 1 \right)}}{\boldsymbol{\dot{w}}} {{\left( x,t \right)}^{T}}{\boldsymbol{\dot{w}}} \left( x,t \right){\rm d}{{A}^{\left( 1 \right)}}{\rm d}x \\ {} & {} & +\frac{1}{2}\int \limits_{{{x}_{e}}}^{{{x}_{e+1}}} \mathop \int \limits_{{{A}^{\left( 2 \right)}}} {{\rho }^{\left( 2 \right)}}{\boldsymbol{\dot{w}}} {{\left( x,t \right)}^{T}}{\boldsymbol{\dot{w}}} \left( x,t \right){\rm d}{{A}^{\left( 2 \right)}}{\rm d}x \\ {} & {} & +I_{0}^{tip}{{x}_{c}}{\boldsymbol{\dot{w}}} {{\left( {{x}_{e+1}},t \right)}^{T}}{\boldsymbol{\dot{ \theta }}} \left( {{x}_{e+1}},t \right) \\ {} & {} & +\frac{1}{2}I_{0}^{tip}{\boldsymbol{\dot{w}}} {{\left( {{x}_{e+1}},t \right)}^{T}}{\boldsymbol{\dot{w}}} \left( {{x}_{e+1}},t \right) \\ {} & {} & +\frac{1}{2}I_{2}^{tip}{\boldsymbol{\dot{ \theta }}} {{\left( {{x}_{e+1}},t \right)}^{T}}{\boldsymbol{\dot{ \theta }}} \left( {{x}_{e+1}},t \right). \\ \end{array}\end{eqnarray} \tag{ 32 }$$

The derivation of equation (32) can be seen in appendix B. The potential energy for the piezoelectric unimorph can be formulated as

$$\begin{eqnarray}\begin{array}{lllllllllllllll} PE & = & \frac{1}{2}\int \limits_{{{x}_{e}}}^{{{x}_{e+1}}} \mathop \int \limits_{{{A}^{\left( 1 \right)}}} {\boldsymbol{S}} _{1}^{{{\left( 1 \right)}^{T}}}{\boldsymbol{T}} _{1}^{\left( 1 \right)}{\rm d}{{A}^{\left( 1 \right)}}{\rm d}x \\ {} & {} & +\frac{1}{2}\int \limits_{{{x}_{e}}}^{{{x}_{e+1}}} \mathop \int \limits_{{{A}^{\left( 2 \right)}}} {\boldsymbol{S}} _{1}^{{{\left( 2 \right)}^{T}}}{\boldsymbol{T}} _{1}^{\left( 2 \right)}{\rm d}{{A}^{\left( 2 \right)}}{\rm d}x. \\ \end{array}\end{eqnarray} \tag{ 33 }$$

The electrical energy term for the piezoelectric layer can be reduced as

$$\begin{eqnarray}&&WE=\frac{1}{2}\int \limits_{{{x}_{e}}}^{{{x}_{e+1}}} \mathop \int \limits_{{{A}^{\left( 2 \right)}}} {\boldsymbol{E}} _{3}^{T}{{{\boldsymbol{D}} }_{3}}{\rm d}{{A}^{\left( 2 \right)}}{\rm d}x.\end{eqnarray} \tag{ 34 }$$

The nonconservative work on the system, due to the input base excitation and electrical charge output, can be written as

$$\begin{eqnarray}\begin{array}{lllllllllllllll} WF & = & -\int \limits_{{{x}_{e}}}^{{{x}_{e+1}}} \mathop \int \limits_{{{A}^{\left( 1 \right)}}} {{\rho }^{\left( 1 \right)}}{\boldsymbol{w}} {{\left( x,t \right)}^{T}}{\rm d}{{A}^{\left( 1 \right)}}{\rm d}x\;{{{{\boldsymbol{\ddot{w}}} }}_{{\boldsymbol{base}} }}\left( t \right) \\ {} & {} & -\int \limits_{{{x}_{e}}}^{{{x}_{e+1}}} \mathop \int \limits_{{{A}^{\left( 2 \right)}}} {{\rho }^{\left( 2 \right)}}{\boldsymbol{w}} {{\left( x,t \right)}^{T}}{\rm d}{{A}^{\left( 2 \right)}}{\rm d}x\;{{{{\boldsymbol{\ddot{w}}} }}_{{\boldsymbol{base}} }}\left( t \right) \\ {} & {} & -I_{0}^{tip}{{x}_{c}}{\boldsymbol{ \theta }} {{\left( {{x}_{e+1}},t \right)}^{T}}{{{{\boldsymbol{\ddot{w}}} }}_{base}}\left( t \right) \\ {} & {} & -I_{0}^{tip}{\boldsymbol{w}} {{\left( {{x}_{e+1}},t \right)}^{T}}{{{{\boldsymbol{\ddot{w}}} }}_{base}}\left( t \right)+{{{\boldsymbol{q}} }^{\left( e \right)}}\left( t \right){{{\boldsymbol{v}} }^{\left( e \right)}}\left( t \right). \\ \end{array}\end{eqnarray} \tag{ 35 }$$

The reduced form of equation (35) due to the input base excitation can be found in appendix C. The expressions of the functional energy forms, including external work, can be modified in terms of the discretized mechanical and electrical shape functions. The extended kinetic energy form can be formulated by substituting equation (21) into equation (32), giving

$$\begin{eqnarray}\begin{array}{lllllllllllllll} KE & = & \frac{1}{2}\int \limits_{{{x}_{e}}}^{{{x}_{e+1}}} \mathop \int \limits_{{{A}^{\left( 1 \right)}}} {{z}^{2}}{{\rho }^{\left( 1 \right)}}{{{{\boldsymbol{\dot{u}}} }}^{e}}{{(t)}^{T}}{{\frac{{\rm d}{{{\boldsymbol{ \Phi }} }^{e}}(x)}{{\rm d}x}}^{T}} \\ {} & {} & \times \frac{{\rm d}{{{\boldsymbol{ \Phi }} }^{e}}(x)}{{\rm d}x}{{{{\boldsymbol{\dot{u}}} }}^{e}}(t){\rm d}{{A}^{\left( 1 \right)}}{\rm d}x \\ {} & {} & +\frac{1}{2}\int \limits_{{{x}_{e}}}^{{{x}_{e+1}}} \mathop \int \limits_{{{A}^{\left( 2 \right)}}} {{z}^{2}}{{\rho }^{\left( 2 \right)}}{{{{\boldsymbol{\dot{u}}} }}^{e}}{{(t)}^{T}}{{\frac{{\rm d}{{{\boldsymbol{ \Phi }} }^{e}}(x)}{{\rm d}x}}^{T}} \\ {} & {} & \times \frac{{\rm d}{{{\boldsymbol{ \Phi }} }^{e}}(x)}{{\rm d}x}{{{{\boldsymbol{\dot{u}}} }}^{e}}(t){\rm d}{{A}^{\left( 2 \right)}}{\rm d}x \\ {} & {} & +\frac{1}{2}\int \limits_{{{x}_{e}}}^{{{x}_{e+1}}} \mathop \int \limits_{{{A}^{\left( 1 \right)}}} {{\rho }^{\left( 1 \right)}}{{{{\boldsymbol{\dot{u}}} }}^{e}}{{(t)}^{T}}{{{\boldsymbol{ \Phi }} }^{e}}{{(x)}^{T}}{{{\boldsymbol{ \Phi }} }^{\left( e \right)}}(x){{{{\boldsymbol{\dot{u}}} }}^{\left( e \right)}}(t){\rm d}{{A}^{\left( 1 \right)}}{\rm d}x \\ {} & {} & +\frac{1}{2}\int \limits_{{{x}_{e}}}^{{{x}_{e+1}}} \mathop \int \limits_{{{A}^{\left( 2 \right)}}} {{\rho }^{\left( 2 \right)}}{{{{\boldsymbol{\dot{u}}} }}^{e}}{{(t)}^{T}}{{{\boldsymbol{ \Phi }} }^{e}}{{(x)}^{T}}{{{\boldsymbol{ \Phi }} }^{e}}(x){{{{\boldsymbol{\dot{u}}} }}^{e}}(t){\rm d}{{A}^{\left( 2 \right)}}{\rm d}x \\ {} & {} & +I_{0}^{tip}{{x}_{c}}{\boldsymbol{\dot{u}}} {{(t)}^{T}}{{{\boldsymbol{ \Phi }} }^{e}}{{\left( {{x}_{e+1}} \right)}^{T}}\frac{d{{{\boldsymbol{ \Phi }} }^{e}}\left( {{x}_{e+1}} \right)}{dx}{\boldsymbol{\dot{u}}} (t) \\ {} & {} & +\frac{1}{2}I_{0}^{tip}{{{{\boldsymbol{\dot{u}}} }}^{e}}{{(t)}^{T}}{{{\boldsymbol{ \Phi }} }^{e}}{{\left( {{x}_{e+1}} \right)}^{T}}{{{\boldsymbol{ \Phi }} }^{e}}\left( {{x}_{e+1}} \right){{{{\boldsymbol{\dot{u}}} }}^{e}}(t) \\ {} & {} & +\frac{1}{2}I_{2}^{tip}{{\frac{d{{{\boldsymbol{ \Phi }} }^{e}}\left( {{x}_{e+1}} \right)}{dx}}^{T}}\;{{{{\boldsymbol{\dot{u}}} }}^{e}}{{(t)}^{T}}\frac{d{{{\boldsymbol{ \Phi }} }^{e}}\left( {{x}_{e+1}} \right)}{dx}{{{{\boldsymbol{\dot{u}}} }}^{e}}(t). \\ \end{array}\end{eqnarray} \tag{ 36 }$$

Corresponding to equations (24) and (29), the extended potential energy form can be reduced as

$$\begin{eqnarray}\begin{array}{lllllllllllllll} PE & = & \frac{1}{2}\int \limits_{{{x}_{e}}}^{{{x}_{e+1}}} \mathop \int \limits_{{{A}^{\left( 1 \right)}}} {{z}^{2}}{{{\boldsymbol{u}} }^{e}}{{(t)}^{T}}{{{\boldsymbol{ \Psi }} }^{e}}{{(x)}^{T}}\bar{c}_{11}^{\left( 1 \right)}{{{\boldsymbol{ \Psi }} }^{e}}(x){{{\boldsymbol{u}} }^{e}}(t){\rm d}{{A}^{\left( 1 \right)}}{\rm d}x \\ {} & {} & +\frac{1}{2}\int \limits_{{{x}_{e}}}^{{{x}_{e+1}}} \mathop \int \limits_{{{A}^{\left( 2 \right)}}} {{z}^{2}}{{{\boldsymbol{u}} }^{e}}{{(t)}^{T}}{{{\boldsymbol{ \Psi }} }^{e}}{{(x)}^{T}}\bar{c}_{11}^{\left( 2 \right)}{{{\boldsymbol{ \Psi }} }^{e}}(x){{{\boldsymbol{u}} }^{e}}(t){\rm d}{{A}^{\left( 2 \right)}}{\rm d}x \\ {} & {} & -\frac{1}{2}\int \limits_{{{x}_{e}}}^{{{x}_{e+1}}} \mathop \int \limits_{{{A}^{\left( 2 \right)}}} z{{e}_{31}}{{{\boldsymbol{u}} }^{e}}{{(t)}^{T}}{{{\boldsymbol{ \Psi }} }^{e}}{{(x)}^{T}}{{{\boldsymbol{ \Omega }} }^{e}}(z){{{\boldsymbol{v}} }^{e}}(t){\rm d}{{A}^{\left( 2 \right)}}{\rm d}x. \\ \end{array}\end{eqnarray} \tag{ 37 }$$

The extended electrical energy expression can be reduced by substituting equations (28) and (30) into equation (34), as

$$\begin{eqnarray}&&\begin{array}{lllllllllllllll} WE=\frac{1}{2}\int \limits_{{{x}_{e}}}^{{{x}_{e+1}}} \mathop \int \limits_{{{A}^{\left( 2 \right)}}} z{{e}_{31}}{{{\boldsymbol{v}} }^{e}}{{(t)}^{T}}{{{\boldsymbol{ \Omega }} }^{e}}{{(z)}^{T}}{{{\boldsymbol{ \Psi }} }^{e}}(x){{{\boldsymbol{u}} }^{e}}(t){\rm d}{{A}^{\left( 2 \right)}}{\rm d}x \\ \ \;+\frac{1}{2}\int \limits_{{{x}_{e}}}^{{{x}_{e+1}}} \mathop \int \limits_{{{A}^{\left( 2 \right)}}} \varepsilon _{33}^{S}{{{\boldsymbol{v}} }^{e}}{{(t)}^{T}}{{{\boldsymbol{ \Omega }} }^{e}}{{(z)}^{T}}{{{\boldsymbol{ \Omega }} }^{e}}(z){{{\boldsymbol{v}} }^{e}}(t){\rm d}{{A}^{\left( 2 \right)}}{\rm d}x. \\ \end{array}\end{eqnarray} \tag{ 38 }$$

The extended virtual work, in terms of the expression of nonconservative energy, can be formulated by substituting equation (21) into (35) to give

$$\begin{eqnarray}\begin{array}{lllllllllllllll} WF & = & -\int \limits_{{{x}_{e}}}^{{{x}_{e+1}}} \mathop \int \limits_{{{A}^{\left( 1 \right)}}} {{\rho }^{\left( 1 \right)}}{{{\boldsymbol{ \Phi }} }^{e}}{{\left( x \right)}^{T}}{{{\boldsymbol{u}} }^{e}}\left( t \right){\rm d}{{A}^{\left( 1 \right)}}{\rm d}x{{{{\boldsymbol{\ddot{w}}} }}_{base}}\left( t \right) \\ {} & {} & -\int \limits_{{{x}_{e}}}^{{{x}_{e+1}}} \mathop \int \limits_{{{A}^{\left( 2 \right)}}} {{\rho }^{\left( 2 \right)}}{{{\boldsymbol{ \Phi }} }^{e}}{{\left( x \right)}^{T}}{{{\boldsymbol{u}} }^{e}}\left( t \right){\rm d}{{A}^{\left( 2 \right)}}{\rm d}x{{{{\boldsymbol{\ddot{w}}} }}_{base}}\left( t \right) \\ {} & {} & -I_{0}^{tip}{{x}_{c}}{{\frac{d{{{\boldsymbol{ \Phi }} }^{e}}\left( {{x}_{e+1}} \right)}{dx}}^{T}}{{{\boldsymbol{u}} }^{e}}\left( t \right){{{{\boldsymbol{\ddot{w}}} }}_{base}}\left( t \right) \\ {} & {} & -I_{0}^{tip}{{{\boldsymbol{ \Phi }} }^{e}}\left( {{x}_{e+1}} \right){{{\boldsymbol{u}} }^{e}}\left( t \right){{{{\boldsymbol{\ddot{w}}} }}_{base}}\left( t \right)+{{{\boldsymbol{q}} }^{e}}\left( t \right){{{\boldsymbol{v}} }^{e}}\left( t \right). \\ \end{array}\end{eqnarray} \tag{ 39 }$$

The expressions given from equations (36)–(39) can be substituted into equation (31) to give two electromechanical dynamic equations. After simplifying, the first electromechanical dynamic equation due to the transverse bending form can be expressed as

$$\begin{eqnarray}&&\begin{array}{lllllllllllllll} \int \limits_{{{x}_{e}}}^{{{x}_{e+1}}} \mathop \int \limits_{{{A}^{\left( 1 \right)}}} {{z}^{2}}{{\frac{{\rm d}{{{\boldsymbol{ \Phi }} }^{e}}\left( x \right)}{{\rm d}x}}^{T}}\frac{{\rm d}{{{\boldsymbol{ \Phi }} }^{e}}\left( x \right)}{{\rm d}x}{{{{\boldsymbol{\ddot{u}}} }}^{e}}\left( t \right){\rm d}{{A}^{\left( 1 \right)}}{\rm d}x \\ \quad +\int \limits_{{{x}_{e}}}^{{{x}_{e+1}}} \mathop \int \limits_{{{A}^{\left( 2 \right)}}} {{z}^{2}}{{\frac{{\rm d}{{{\boldsymbol{ \Phi }} }^{e}}\left( x \right)}{{\rm d}x}}^{T}}\frac{{\rm d}{{{\boldsymbol{ \Phi }} }^{e}}\left( x \right)}{{\rm d}x}{{{{\boldsymbol{\ddot{u}}} }}^{e}}\left( t \right){\rm d}{{A}^{\left( 2 \right)}}{\rm d}x \\ \quad +\int \limits_{{{x}_{e}}}^{{{x}_{e+1}}} \mathop \int \limits_{{{A}^{\left( 1 \right)}}} {{\rho }^{\left( 1 \right)}}{{{\boldsymbol{ \Phi }} }^{e}}{{\left( x \right)}^{T}}{{{\boldsymbol{ \Phi }} }^{e}}\left( x \right){{{{\boldsymbol{\ddot{u}}} }}^{e}}\left( t \right){\rm d}{{A}^{\left( 1 \right)}}{\rm d}x \\ \quad +\int \limits_{{{x}_{e}}}^{{{x}_{e+1}}} \mathop \int \limits_{{{A}^{\left( 2 \right)}}} {{\rho }^{\left( 2 \right)}}{{{\boldsymbol{ \Phi }} }^{e}}{{\left( x \right)}^{T}}{{{\boldsymbol{ \Phi }} }^{\left( e \right)}}\left( x \right){{{{\boldsymbol{\ddot{u}}} }}^{\left( e \right)}}\left( t \right){\rm d}{{A}^{\left( 2 \right)}}{\rm d}x \\ \quad +2I_{0}^{tip}{{x}_{c}}{{{\boldsymbol{ \Phi }} }^{e}}{{\left( {{x}_{e+1}} \right)}^{T}}\frac{{\rm d}{{{\boldsymbol{ \Phi }} }^{e}}\left( {{x}_{e+1}} \right)}{{\rm d}x}{\boldsymbol{\ddot{u}}} \left( t \right) \\ \quad +I_{0}^{tip}{{{\boldsymbol{ \Phi }} }^{e}}{{\left( {{x}_{e+1}} \right)}^{T}}{{{\boldsymbol{ \Phi }} }^{e}}\left( {{x}_{e+1}} \right){{{{\boldsymbol{\ddot{u}}} }}^{e}}\left( t \right) \\ \quad +I_{2}^{tip}{{\frac{{\rm d}{{{\boldsymbol{ \Phi }} }^{e}}\left( {{x}_{e+1}} \right)}{{\rm d}x}}^{T}}\frac{{\rm d}{{{\boldsymbol{ \Phi }} }^{e}}\left( {{x}_{e+1}} \right)}{{\rm d}x}{{{{\boldsymbol{\ddot{u}}} }}^{e}}\left( t \right) \\ \quad +\int \limits_{{{x}_{e}}}^{{{x}_{e+1}}} \mathop \int \limits_{{{A}^{\left( 1 \right)}}} {{z}^{2}}\bar{c}_{11}^{\left( 1 \right)}{{{\boldsymbol{ \Psi }} }^{e}}\left( x \right){{{\boldsymbol{ \Psi }} }^{e}}{{\left( x \right)}^{T}}{{{\boldsymbol{u}} }^{e}}\left( t \right){\rm d}{{A}^{\left( 1 \right)}}{\rm d}x \\ \quad +\int \limits_{{{x}_{e}}}^{{{x}_{e+1}}} \mathop \int \limits_{{{A}^{\left( 2 \right)}}} {{z}^{2}}\bar{c}_{11}^{\left( 2 \right)}{{{\boldsymbol{ \Psi }} }^{e}}\left( x \right){{{\boldsymbol{ \Psi }} }^{e}}{{\left( x \right)}^{T}}{{{\boldsymbol{u}} }^{e}}\left( t \right){\rm d}{{A}^{\left( 2 \right)}}{\rm d}x \\ \quad -\frac{1}{2}\int \limits_{{{x}_{e}}}^{{{x}_{e+1}}} \mathop \int \limits_{{{A}^{\left( 2 \right)}}} z{{e}_{31}}{{{\boldsymbol{ \Psi }} }^{e}}{{\left( x \right)}^{T}}{{{\boldsymbol{ \Omega }} }^{e}}\left( z \right){{{\boldsymbol{v}} }^{e}}\left( t \right){\rm d}{{A}^{\left( 2 \right)}}{\rm d}x \\ \quad -\frac{1}{2}\int \limits_{{{x}_{e}}}^{{{x}_{e+1}}} \mathop \int \limits_{{{A}^{\left( 2 \right)}}} z{{e}_{31}}{{{\boldsymbol{ \Psi }} }^{e}}{{\left( x \right)}^{T}}{{{\boldsymbol{ \Omega }} }^{e}}\left( z \right){{{\boldsymbol{v}} }^{e}}\left( t \right){\rm d}{{A}^{\left( 2 \right)}}{\rm d}x \\ =-\int \limits_{{{x}_{e}}}^{{{x}_{e+1}}} \mathop \int \limits_{{{A}^{\left( 1 \right)}}} {{\rho }^{\left( 1 \right)}}{{{\boldsymbol{ \Phi }} }^{e}}{{\left( x \right)}^{T}}{\rm d}{{A}^{\left( 1 \right)}}{\rm d}x{{{{\boldsymbol{\ddot{w}}} }}_{base}}\left( t \right) \\ \quad -\int \limits_{{{x}_{e}}}^{{{x}_{e+1}}} \mathop \int \limits_{{{A}^{\left( 2 \right)}}} {{\rho }^{\left( 2 \right)}}{{{\boldsymbol{ \Phi }} }^{e}}{{\left( x \right)}^{T}}{\rm d}{{A}^{\left( 2 \right)}}{\rm d}x{{{{\boldsymbol{\ddot{w}}} }}_{base}}\left( t \right) \\ \qquad \qquad \qquad \qquad \quad -I_{0}^{tip}{{x}_{c}}{{\frac{{\rm d}{{{\boldsymbol{ \Phi }} }^{e}}({{L}_{e}})}{{\rm d}x}}^{T}}{{{{\boldsymbol{\ddot{w}}} }}_{base}}\left( t \right)-I_{0}^{tip}{{{\boldsymbol{ \Phi }} }^{e}}{{\left( {{L}_{e}} \right)}^{T}}{{{{\boldsymbol{\ddot{w}}} }}_{base}}\left( t \right). \\ \end{array}\end{eqnarray} \tag{ 40 }$$

The second electromechanical dynamic equation due to the electrical form can be expressed as

$$\begin{eqnarray}&&\begin{array}{lllllllllllllll} -\frac{1}{2}\int \limits_{{{x}_{e}}}^{{{x}_{e+1}}} \mathop \int \limits_{{{A}^{\left( 2 \right)}}} z{{e}_{31}}{{{\boldsymbol{ \Omega }} }^{e}}{{\left( z \right)}^{T}}{{{\boldsymbol{ \Psi }} }^{e}}\left( x \right){{{\boldsymbol{u}} }^{e}}\left( t \right){\rm d}{{A}^{\left( 2 \right)}}{\rm d}x \\ \quad -\frac{1}{2}\int \limits_{{{x}_{e}}}^{{{x}_{e+1}}} \mathop \int \limits_{{{A}^{\left( 2 \right)}}} z{{e}_{31}}{{{\boldsymbol{ \Omega }} }^{e}}{{\left( z \right)}^{T}}{{{\boldsymbol{ \Psi }} }^{e}}\left( x \right){{{\boldsymbol{u}} }^{e}}\left( t \right){\rm d}{{A}^{\left( 2 \right)}}{\rm d}x \\ \quad -\int \limits_{{{x}_{e}}}^{{{x}_{e+1}}} \mathop \int \limits_{{{A}^{\left( 2 \right)}}} \varepsilon _{33}^{S}{\boldsymbol{ \Omega }} {{\left( z \right)}^{T}}{\boldsymbol{ \Omega }} \left( z \right){\boldsymbol{v}} \left( t \right){\rm d}{{A}^{\left( 2 \right)}}{\rm d}x-{{{\boldsymbol{q}} }^{e}}\left( t \right)=0. \\ \end{array}\end{eqnarray} \tag{ 41 }$$

Equation (41) can be modified by differentiating with respect to time, to give

$$\begin{eqnarray}&&\begin{array}{lllllllllllllll} -\frac{1}{2}\int \limits_{{{x}_{e}}}^{{{x}_{e+1}}} \mathop \int \limits_{{{A}^{\left( 2 \right)}}} z{\boldsymbol{ \Omega }} {{\left( z \right)}^{T}}{{e}_{31}}{\boldsymbol{ \Psi }} \left( x \right){\boldsymbol{\dot{u}}} \left( t \right){\rm d}{{A}^{\left( 2 \right)}}{\rm d}x \\ \quad -\frac{1}{2}\int \limits_{{{x}_{e}}}^{{{x}_{e+1}}} \mathop \int \limits_{{{A}^{\left( 2 \right)}}} z{\boldsymbol{ \Omega }} {{\left( z \right)}^{T}}{{e}_{31}}{\boldsymbol{ \Psi }} \left( x \right){\boldsymbol{\dot{u}}} \left( t \right){\rm d}{{A}^{\left( 2 \right)}}{\rm d}x \\ \quad -\int \limits_{{{x}_{e}}}^{{{x}_{e+1}}} \mathop \int \limits_{{{A}^{\left( 2 \right)}}} \varepsilon _{33}^{S}{\boldsymbol{ \Omega }} {{\left( z \right)}^{T}}{\boldsymbol{ \Omega }} \left( z \right){\boldsymbol{\dot{v}}} \left( t \right){\rm d}{{A}^{\left( 2 \right)}}{\rm d}x{\mkern 1mu} -{\boldsymbol{i}} _{p}^{e}\left( t \right)=0. \\ \end{array}\end{eqnarray} \tag{ 42 }$$

4.3. Element matrices of electromechanical dynamic equations

Constitutive electromechanical nonhomogeneous differential dynamic equations, in terms of equations (40) and (42), can be arranged into matrix form by including Rayleigh damping to give

$$\begin{eqnarray}\begin{array}{lllllllllllllll} \left[ \begin{array}{ccccccccccccccc} {{{\boldsymbol{M}} }^{e}} & 0 \\ 0 & 0 \\ \end{array} \right]\left\{ \begin{array}{ccccccccccccccc} {{{{\boldsymbol{\ddot{u}}} }}^{e}}\left( t \right) \\ {{{{\boldsymbol{\ddot{v}}} }}^{e}}\left( t \right) \\ \end{array} \right\}+\left[ \begin{array}{ccccccccccccccc} {{{\boldsymbol{C}} }^{e}} & 0 \\ {\boldsymbol{P}} _{\theta }^{eT} & {\boldsymbol{P}} _{D}^{e} \\ \end{array} \right]\left\{ \begin{array}{ccccccccccccccc} {{{{\boldsymbol{\dot{u}}} }}^{e}}\left( t \right) \\ {{{{\boldsymbol{\dot{v}}} }}^{e}}\left( t \right) \\ \end{array} \right\} \\ \quad +\left[ \begin{array}{ccccccccccccccc} {{{\boldsymbol{K}} }^{e}} & {\boldsymbol{P}} _{\theta }^{e} \\ 0 & 0 \\ \end{array} \right]\left\{ \begin{array}{ccccccccccccccc} {{{\boldsymbol{u}} }^{e}}\left( t \right) \\ {{{\boldsymbol{v}} }^{e}}\left( t \right) \\ \end{array} \right\}=\left\{ \begin{array}{ccccccccccccccc} {{{\boldsymbol{F}} }^{e}} \\ {\boldsymbol{i}} _{p}^{e} \\ \end{array} \right\}, \\ \end{array}\end{eqnarray} \tag{ 43 }$$

where

$$\begin{eqnarray}\begin{array}{lllllllllllllll} {{{\boldsymbol{M}} }^{e}} & = & \mathop \sum \limits_{nl=1} \left( \int \limits_{{{x}_{e}}}^{{{x}_{e+1}}} \mathop \int \limits_{{{A}^{\left( nl \right)}}} {{z}^{2}}{{\rho }^{\left( nl \right)}}{{\frac{{\rm d}{{{\boldsymbol{ \Phi }} }^{e}}\left( x \right)}{{\rm d}x}}^{T}}\frac{{\rm d}{{{\boldsymbol{ \Phi }} }^{e}}\left( x \right)}{{\rm d}x}{\rm d}{{A}^{\left( nl \right)}}{\rm d}x \right) \\ {} & {} & +\mathop \sum \limits_{nl=1} \left( \int \limits_{{{x}_{e}}}^{{{x}_{e+1}}} \mathop \int \limits_{{{A}^{\left( nl \right)}}} {{\rho }^{\left( nl \right)}}{{{\boldsymbol{ \Phi }} }^{e}}{{\left( x \right)}^{T}}{{{\boldsymbol{ \Phi }} }^{e}}\left( x \right){\rm d}{{A}^{\left( nl \right)}}{\rm d}x \right) \\ {} & {} & +{\mkern 1mu} 2I_{0}^{tip}{{x}_{c}}{{{\boldsymbol{ \Phi }} }^{e}}{{\left( {{x}_{e+1}} \right)}^{T}}\frac{{\rm d}{{{\boldsymbol{ \Phi }} }^{e}}\left( {{x}_{e+1}} \right)}{{\rm d}x} \\ {} & {} & +{\mkern 1mu} I_{0}^{tip}{{{\boldsymbol{ \Phi }} }^{e}}{{\left( {{x}_{e+1}} \right)}^{T}}{{{\boldsymbol{ \Phi }} }^{e}}\left( {{x}_{e+1}} \right) \\ {} & {} & +{\mkern 1mu} I_{2}^{tip}{{\frac{{\rm d}{{{\boldsymbol{ \Phi }} }^{e}}\left( {{x}_{e+1}} \right)}{{\rm d}x}}^{T}}\frac{{\rm d}{{{\boldsymbol{ \Phi }} }^{e}}\left( {{x}_{e+1}} \right)}{{\rm d}x} \\ {{{\boldsymbol{K}} }^{e}} & = & \mathop \sum \limits_{nl=1} \left( \int \limits_{{{x}_{e}}}^{{{x}_{e+1}}} \mathop \int \limits_{{{A}^{\left( nl \right)}}} {{z}^{2}}\bar{c}_{11}^{\left( nl \right)}{{{\boldsymbol{ \Psi }} }^{e}}{{\left( x \right)}^{T}}{{{\boldsymbol{ \Psi }} }^{e}}\left( x \right){\rm d}{{A}^{\left( nl \right)}}{\rm d}x \right), \\ {{{\boldsymbol{C}} }^{e}} & = & \alpha {{{\boldsymbol{M}} }^{e}}+\beta {{{\boldsymbol{K}} }^{e}} \\ {\boldsymbol{P}} _{\theta }^{e} & = & -\mathop \sum \limits_{np=1} \left( \int \limits_{{{x}_{e}}}^{{{x}_{e+1}}} \mathop \int \limits_{{{A}^{\left( np \right)}}} ze_{31}^{\left( np \right)}{{{\boldsymbol{ \Psi }} }^{e}}{{\left( x \right)}^{T}}{{{\boldsymbol{ \Omega }} }^{e}}\left( z \right){\rm d}{{A}^{\left( np \right)}}{\rm d}x \right), \\ {\boldsymbol{P}} _{\theta }^{{{e}^{T}}} & = & -\mathop \sum \limits_{np=1} \left( \int \limits_{{{x}_{e}}}^{{{x}_{e+1}}} \mathop \int \limits_{{{A}^{\left( np \right)}}} ze_{31}^{\left( np \right)}{{{\boldsymbol{ \Omega }} }^{e}}{{\left( z \right)}^{T}}{{{\boldsymbol{ \Psi }} }^{e}}\left( x \right){\rm d}{{A}^{\left( np \right)}}{\rm d}x \right) \\ {\boldsymbol{P}} _{D}^{e} & = & -\mathop \sum \limits_{np=1} \left( \int \limits_{{{x}_{e}}}^{{{x}_{e+1}}} \mathop \int \limits_{{{A}^{\left( np \right)}}} \varepsilon _{33}^{\left( np \right)}{{{\boldsymbol{ \Omega }} }^{e}}{{\left( z \right)}^{T}}{{{\boldsymbol{ \Omega }} }^{e}}\left( z \right){\rm d}{{A}^{\left( np \right)}}{\rm d}x \right), \\ {{{\boldsymbol{F}} }^{e}} & = & -{{{\boldsymbol{Q}} }^{e}}{{{{\boldsymbol{\ddot{w}}} }}_{base}}\left( t \right) \\ {{{\boldsymbol{Q}} }^{\left( e \right)}} & = & \mathop \sum \limits_{nl=1} \left( \int \limits_{{{x}_{e}}}^{{{x}_{e+1}}} \mathop \int \limits_{{{A}^{\left( nl \right)}}} {{\rho }^{\left( nl \right)}}{{{\boldsymbol{ \Phi }} }^{e}}{{\left( x \right)}^{T}}{\rm d}{{A}^{\left( nl \right)}}{\rm d}x \right) \\ {} & {} & +I_{0}^{tip}{{x}_{c}}{{\frac{{\rm d}{{{\boldsymbol{ \Phi }} }^{e}}\left( {{x}_{e+1}} \right)}{{\rm d}x}}^{T}}+I_{0}^{tip}{{{\boldsymbol{ \Phi }} }^{e}}{{\left( {{x}_{e+1}} \right)}^{T}}. \\ \end{array}\end{eqnarray} \tag{ 44 }$$

It is important to note here that equation (43) consists of mass, stiffness, electromechanical coupling and piezoelectric capacitance matrices and dynamic force vectors. These parameters are modelled at the local electromechanical element, which is located between the end of the unimorph beam length and the tip mass. The contribution of the offset tip mass can be seen in the mass matrices and dynamic force vectors, which most of the existing literature has ignored. As the unimorph-based Euler-Bernoulli beam assumption is considered, the first term from the mass matrix (rotary inertia of the unimorph) can be discarded due to the contribution of the Rayleigh beam assumption. Also, the contribution of input base excitation to the cantilever unimorph beam mathematically affects the distributed input dynamic forces into all the elements. Note that when other elements of the unimorph beam are considered in the numerical simulation, the parameter for the tip mass should be discarded from the mass matrices and input dynamic force vectors.

4.4. Normalized global element matrices of electromechanical dynamic equations

Constitutive electromechanical finite element dynamic equations in global matrix form reduced from the global transformation of local elements can be formulated as

$$\begin{eqnarray}&&\begin{array}{lllllllllllllll} \mathop{{\boldsymbol{M}} }\limits_{\left( {{n}_{m}}\times {{n}_{m}} \right)}\mathop{{\boldsymbol{\ddot{u}}} \left( t \right)}\limits_{\left( {{n}_{m}}\times 1 \right)}+\mathop{{\boldsymbol{C}} }\limits_{\left( {{n}_{m}}\times {{n}_{m}} \right)}\mathop{{\boldsymbol{\dot{u}}} \left( t \right)}\limits_{\left( {{n}_{m}}\times 1 \right)}+\mathop{{\boldsymbol{K}} }\limits_{\left( {{n}_{m}}\times {{n}_{m}} \right)}\mathop{{\boldsymbol{u}} \left( t \right)}\limits_{\left( {{n}_{m}}\times 1 \right)}+\mathop{{{{\boldsymbol{P}} }_{\theta }}}\limits_{\left( {{n}_{m}}\times {{n}_{e}} \right)}\mathop{{\boldsymbol{v}} \left( t \right)}\limits_{\left( {{n}_{e}}\times 1 \right)} \\ \qquad =\mathop{{\boldsymbol{F}} }\limits_{\left( {{n}_{m}}\times 1 \right)}\mathop{{\boldsymbol{P}} _{\theta }^{T}}\limits_{\left( {{n}_{e}}\times {{n}_{m}} \right)}\mathop{{\boldsymbol{\dot{u}}} \left( t \right)}\limits_{\left( {{n}_{m}}\times 1 \right)}+\mathop{{{{\boldsymbol{P}} }_{D}}}\limits_{\left( {{n}_{e}}\times {{n}_{e}} \right)}\mathop{{\boldsymbol{\dot{v}}} \left( t \right)}\limits_{\left( {{n}_{e}}\times 1 \right)}=\mathop{{{{\boldsymbol{i}} }_{p}}\left( t \right)}\limits_{\left( {{n}_{e}}\times 1 \right)}. \\ \end{array}\end{eqnarray} \tag{ 45 }$$

Note that the script terms ${{n}_{m}}$ and ${{n}_{e}}$ showing at each parameter in equation (45) indicate the global matrices corresponding to global displacement vectors at nodes (nodal degrees of freedom) and global electrical voltage vectors at elements (electrical degrees of freedom), respectively. For assembling the global matrices, the connectivity matrices whose discretized elements share the same node should meet the requirements of compatibility. As the unimorph beam is a cantilevered structure under input base excitation, the fixed boundary condition at the base structure must be employed in the global matrices. For this case, the derivations for the normalized electromechanical dynamic equations and FRFs can be solved and given in the next stage. The use of the MATLAB computer program not only provides the best method for solving this case, but also gives an efficient computation technique.

The solution form of time-dependent displacement can be stated in terms of the normalized modal vector and the time-dependent displacement generalized coordinate as

$$\begin{eqnarray}\begin{array}{lllllllllllllll} {\boldsymbol{u}} \left( t \right) & = & {{{\boldsymbol{ \varphi }} }_{1}}{{{\boldsymbol{a}} }_{1}}\left( t \right)+{{{\boldsymbol{ \varphi }} }_{2}}{{{\boldsymbol{a}} }_{2}}\left( t \right)+... \\ {} & {} & +{{{\boldsymbol{ \varphi }} }_{m-1}}{{{\boldsymbol{a}} }_{m-1}}\left( t \right)+{{{\boldsymbol{ \varphi }} }_{m}}{{{\boldsymbol{a}} }_{m}}\left( t \right)={\boldsymbol{ \varphi a}} \left( t \right), \\ \end{array}\end{eqnarray} \tag{ 46 }$$

where the normalized modal matrix can be formulated as,

$$\begin{eqnarray}&&{\boldsymbol{ \varphi }} ={\boldsymbol{U}} {{\left( {{{\boldsymbol{U}} }^{T}}{\boldsymbol{MU}} \right)}^{-\,\frac{1}{2}}}.\end{eqnarray} \tag{ 47 }$$

It should be noted that parameters ${\boldsymbol{ \varphi }}$ and ${\boldsymbol{U}}$ represent a set of the normalized modal matrix and generalized eigenvector, respectively. To obtain the eigenvector and the eigenvalue solution, the undamped mechanical dynamic equation from the stiffness and mass matrices should be nonsingular, since the eigenvalue equation, $\left( K-{{\omega }^{2}}M \right)U=0$, can be used to find the nontrivial eigenvector. Note that the normalized modal matrix is the matrix whose columns consist of the normalized eigenvectors [56]. The elemental strain field of the unimorph beam given in equation (24) can be formulated in terms of equations (24) and (46) to give

$$\begin{eqnarray}\begin{array}{lllllllllllllll} {\boldsymbol{S}} \left( x,t \right) & = & -z{{{\boldsymbol{ \Psi }} }^{e}}\left( x \right){{{\boldsymbol{u}} }^{e}}\left( t \right)\mathop{=-z{{{\boldsymbol{ \Psi }} }^{e}}\left( x \right)}\limits_{\left( 1\times 4 \right)}\mathop{{{{\boldsymbol{ \varphi }} }^{e}}}\limits_{\left( 4\times r \right)}\mathop{{\boldsymbol{a}} \left( t \right)}\limits_{\left( r\times 1 \right)} \\ {} & = & \mathop{-z{{{\boldsymbol{ \mu }} }^{e}}\left( x \right)}\limits_{\left( 1\times r \right)}\mathop{{\boldsymbol{a}} \left( t \right)}\limits_{\left( r\times 1 \right)}. \\ \end{array}\end{eqnarray} \tag{ 48 }$$

It should be noted that the elemental strain node parameter, ${{\mu }^{e}}\left( x \right)$, based on the normalized eigenvector, can be determined at certain elements and nodes of the system over the interval, ${{x}_{e}}\leqslant x\leqslant {{x}_{e+1}}$. Note that collective data points for strain nodes along the unimorph beam can only be identified for each discretized element using the differential strain–displacement function, ${{\Psi }^{e}}\left( x \right)$, and elemental normalized eigenvector, ${\boldsymbol{ \varphi }}$. The normalized eigenvector at certain elements can only be determined by separating the cell vector array from the normalized modal matrix in terms of the elemental degree of freedom and the multimode system.

Equation (47) can be used to diagonalize the mass, stiffness and damping matrices to simplify the computation technique. Substituting equation (46) into equation (45) and premultiplying the result by ${{\varphi }^{T}}$ gives

$$\begin{eqnarray}&&\begin{array}{lllllllllllllll} {{{\boldsymbol{ \varphi }} }^{T}}{\bf M}{\boldsymbol{ \varphi \ddot{a}}} \left( t \right)+{{{\boldsymbol{ \varphi }} }^{T}}{\rm C}{\boldsymbol{ \varphi \dot{a}}} \left( t \right)+{{{\boldsymbol{ \varphi }} }^{T}}K{\boldsymbol{ \varphi a}} \left( t \right)+{{{\boldsymbol{ \varphi }} }^{T}}{\mkern 1mu} {{{\boldsymbol{P}} }_{\theta }}{\boldsymbol{v}} \left( t \right)={{{\boldsymbol{ \varphi }} }^{T}}{\boldsymbol{F}} , \\ {\boldsymbol{P}} _{\theta }^{T}{\boldsymbol{ \varphi \dot{a}}} \left( t \right)+{{{\boldsymbol{P}} }_{D}}{\boldsymbol{\dot{v}}} \left( t \right)={{{\boldsymbol{i}} }_{p}}\left( t \right), \\ \end{array}\end{eqnarray} \tag{ 49 }$$

or simplifying equation (49) becomes

$$\begin{eqnarray}&&\begin{array}{lllllllllllllll} {\boldsymbol{\ddot{a}}} \left( t \right)+2{\boldsymbol{ \zeta \omega \dot{a}}} \left( t \right)+{{{\boldsymbol{ \omega }} }^{2}}{\boldsymbol{a}} \left( t \right)+{{{{\boldsymbol{\hat{P}}} }}_{\theta }}{\boldsymbol{v}} \left( t \right)=-\hat{Q}{{{{\boldsymbol{\ddot{w}}} }}_{base}}\left( t \right) \\ \quad {\boldsymbol{\hat{P}}} _{\theta }^{T}{\boldsymbol{\dot{a}}} \left( t \right)+{{{\boldsymbol{P}} }_{D}}{\boldsymbol{\dot{v}}} \left( t \right)={{{\boldsymbol{i}} }_{p}}\left( t \right). \\ \end{array}\end{eqnarray} \tag{ 50 }$$

Parameter $\hat{Q}{{\ddot{w}}_{base}}\left( t \right)$ from the first line of equation (50) indicates the nodal dynamic forces of the discretized unimorph harvester model due to the input base excitation. Other normalized parameters from equation (50) can be stated as

$$\begin{eqnarray}\begin{array}{lllllllllllllll} {{{\boldsymbol{ \varphi }} }^{T}}M{\boldsymbol{ \varphi }} & = & I,\;{{{\boldsymbol{ \varphi }} }^{T}}K{\boldsymbol{ \varphi }} ={{{\boldsymbol{ \omega }} }^{2}}, \\ {{{\boldsymbol{ \varphi }} }^{T}}{\rm C}{\boldsymbol{ \varphi }} & = & {\bf \alpha }\left( {{{\boldsymbol{ \varphi }} }^{T}}{\rm M}{\boldsymbol{ \varphi }} \right)+\beta \left( {{{\boldsymbol{ \varphi }} }^{T}}K{\boldsymbol{ \varphi }} \right)=\alpha {\rm I}+\beta {{{\boldsymbol{ \omega }} }^{2}}=2\zeta {\boldsymbol{ \omega }} \\ {} & {} & {{{{\boldsymbol{\hat{P}}} }}_{\theta }}={{{\boldsymbol{ \varphi }} }^{T}}{{{\boldsymbol{P}} }_{\theta }},\quad {\boldsymbol{\hat{P}}} _{{\boldsymbol{ \theta }} }^{{\rm T}}={\boldsymbol{P}} _{{\boldsymbol{ \theta }} }^{{\rm T}}\varphi ,\quad {\boldsymbol{\hat{Q}}} ={{{\boldsymbol{ \varphi }} }^{T}}{\boldsymbol{Q}} . \\ \end{array}\end{eqnarray} \tag{ 51 }$$

It should be noted that the first two lines of equation (51) represent the orthornormality property of mechanical dynamic equations, the results of which indicate diagonal matrices.

4.5. Normalized global scalar form of the electromechanical dynamic equations

In simplification, the normalized global element matrices of electromechanical transverse equations from equation (50) can be reduced into scalar forms in order to formulate the multimode frequency analysis of the distributed piezoelectric component, as shown in figure 2(a). The electromechanical equations consist of the coupled field normalized differential equations for each elemental global coordinate component. In this case, the first form of the discretized electromechanical piezoelectric dynamic equation can be formulated for the multidegree of freedom (multimode) system, $r=1,2,3,\ldots ,NDOF$, in terms of the number of normalized piezoelectric elements, $s=1,2,3,\ldots ,NELP$, as

$$\begin{eqnarray}&&\begin{array}{lllllllllllllll} {{\ddot{a}}_{r}}\left( t \right)+2{{\zeta }_{r}}{{\omega }_{r}}{{{\dot{a}}}_{{\rm r}}}\left( {\rm t} \right)+\omega _{r}^{2}{{a}_{r}}\left( t \right)+\sum \limits_{s=1}^{NELP} {{{\hat{P}}}_{rs}}{{v}_{s}}\left( t \right) \\ \quad =-{{{\hat{Q}}}_{r}}{{{\ddot{w}}}_{base}}\left( t \right), \\ r=1,2,\ldots ,NDOF. \\ \end{array}\end{eqnarray} \tag{ 52 }$$

The second form of the discretized electromechanical piezoelectric dynamic equation can be formulated as

$$\begin{eqnarray}&&\begin{array}{lllllllllllllll} \sum \limits_{r=1}^{NDOF} {{{\hat{P}}}_{sr}}{{{\dot{a}}}_{r}}\left( t \right)+{{P}_{Ds}}{{{\dot{v}}}_{s}}\left( t \right)={{i}_{Ps}}\left( t \right), \\ \quad s=1,2,3,\ldots ,NELP. \\ \end{array}\end{eqnarray} \tag{ 53 }$$

Note that since the piezoelectric surface with the electrode layers is distributed evenly on the substructure layer, the coupled equation for the global coordinate system should be based on the electromechanical discretization, as shown in figure 2(b). The internal parallel connection, in terms of the Kirchhoff's voltage law and Kirchhoff's current law, can be used to formulate the electrical discretized elements as

$$\begin{eqnarray}&&\begin{array}{lllllllllllllll} {{v}_{1}}\left( t \right)={{v}_{2}}\left( t \right)=....={{v}_{s}}\left( t \right)=v\left( t \right) \\ {{i}_{P1}}\left( t \right)+{{i}_{P2}}\left( t \right)+....+{{i}_{Ps}}\left( t \right)={{i}_{P}}\left( t \right). \\ \end{array}\end{eqnarray} \tag{ 54 }$$

For voltage output, the external load resistance can be used as an external circuit to connect all electrical discretized piezoelectric elements in a parallel connection, as

$$\begin{eqnarray}&&v\left( t \right)={{i}_{P}}\left( t \right){{R}_{load}}.\end{eqnarray} \tag{ 55 }$$

Multimode FRFs of the distributed piezoelectric unimorph can be formulated in several ways. First, one can modify the first term of equation (53), which algebraically corresponds to the number of the normalized piezoelectric elements. Second, one can employ equations (54) and (55) in the equations obtained from the first step. Third, one can apply equation (54) in equation (52). Fourth, the results obtained from the second and third steps can be algebraically solved using Laplace transforms, giving the result in matrix form. The superposition of equations for voltage multimode FRFs can be formulated after simplifying, as

$$\begin{eqnarray}&&\begin{array}{lllllllllllllll} \frac{v\left( j\omega \right)}{-{{\omega }^{2}}{{w}_{base}}{{e}^{j\omega t}}} \\ \quad =\frac{\sum \limits_{r=1}^{NDOF} \frac{{{{\hat{Q}}}_{r}}j\omega \sum \limits_{s=1}^{NELP} {{{\hat{P}}}_{sr}}}{\omega _{r}^{2}-{{\omega }^{2}}+j2{{\zeta }_{r}}{{\omega }_{r}}\omega }}{\sum \limits_{s=1}^{NELP} j\omega {{P}_{Ds}}-\frac{1}{{{R}_{load}}}-\sum \limits_{r=1}^{NDOF} \frac{j\omega \sum \limits_{s=1}^{NELP} {{{\hat{P}}}_{sr}}\sum \limits_{s=1}^{NELP} {{{\hat{P}}}_{rs}}}{\omega _{r}^{2}-{{\omega }^{2}}+j2{{\zeta }_{r}}{{\omega }_{r}}\omega }}. \\ \end{array}\end{eqnarray} \tag{ 56 }$$

The multimode FRF of the electric current output related to the input base transverse acceleration can be stated as

$$\begin{eqnarray}&&\begin{array}{lllllllllllllll} \frac{{{i}_{p}}\left( j\omega \right)}{-{{\omega }^{2}}{{w}_{base}}{{e}^{j\omega t}}} \\ \quad =\frac{\frac{1}{{{R}_{load}}}\sum \limits_{r=1}^{NDOF} \frac{{{{\hat{Q}}}_{r}}j\omega \sum \limits_{s=1}^{NELP} {{{\hat{P}}}_{sr}}}{\omega _{r}^{2}-{{\omega }^{2}}+j2{{\zeta }_{r}}{{\omega }_{r}}\omega }}{\sum \limits_{s=1}^{NELP} j\omega {{P}_{Ds}}-\frac{1}{{{R}_{load}}}-\sum \limits_{r=1}^{NDOF} \frac{j\omega \sum \limits_{s=1}^{NELP} {{{\hat{P}}}_{sr}}\sum \limits_{s=1}^{NELP} {{{\hat{P}}}_{rs}}}{\omega _{r}^{2}-{{\omega }^{2}}+j2{{\zeta }_{r}}{{\omega }_{r}}\omega }.} \\ \end{array}\end{eqnarray} \tag{ 57 }$$

The power harvesting multimode FRF, related to the input transverse acceleration, can be formulated as

$$\begin{eqnarray}&&\begin{array}{lllllllllllllll} \frac{{{P}_{p}}\left( j\omega \right)}{{{\left( -{{\omega }^{2}}{{w}_{base}}{{e}^{j\omega t}} \right)}^{2}}} \\ ={{\left[ \frac{\begin{array}{lllllllllllllll} \left[ \frac{1}{\sqrt{{{R}_{load}}}}\sum \limits_{r=1}^{NDOF} \frac{{{{\hat{Q}}}_{r}}j\omega \sum \limits_{s=1}^{NELP} {{{\hat{P}}}_{sr}}}{\omega _{r}^{2}-{{\omega }^{2}}+j2{{\zeta }_{r}}{{\omega }_{r}}\omega } \right] \\ \end{array}}{\begin{array}{lllllllllllllll} \left[ \sum \limits_{s=1}^{NELP} j\omega {{P}_{Ds}}-\frac{1}{{{R}_{load}}}-\sum \limits_{r=1}^{NDOF} \frac{j\omega \sum \limits_{s=1}^{NELP} {{{\hat{P}}}_{sr}}\sum \limits_{s=1}^{NELP} {{{\hat{P}}}_{rs}}}{\omega _{r}^{2}-{{\omega }^{2}}+j2{{\zeta }_{r}}{{\omega }_{r}}\omega } \right] \\ \end{array}} \right]}^{2}}. \\ \end{array}\end{eqnarray} \tag{ 58 }$$

To formulate the optimal multimode FRF of power harvesting, the reduced optimal load resistance needs to be formulated first. For this case, equation (58) can be differentiated with respect to load resistance, and the differentiable power function can be set to zero to give the optimal load resistance as

$$\begin{eqnarray}&&R_{load}^{opt}=\frac{\sqrt{X{{\left( \omega \right)}^{2}}+Y{{\left( \omega \right)}^{2}}}}{X{{\left( \omega \right)}^{2}}+Y{{\left( \omega \right)}^{2}}},\end{eqnarray} \tag{ 59 }$$

where

$$\begin{eqnarray}&&\begin{array}{lllllllllllllll} X\left( \omega \right)=\sum \limits_{{\rm s}=1}^{NELP} \omega {{P}_{{\rm Ds}}}-\sum \limits_{r=1}^{NDOF} \frac{\omega \sum \limits_{s=1}^{NELP} {{{\hat{P}}}_{{\rm sr}}}\sum \limits_{s=1}^{NELP} {{{\hat{P}}}_{rs}}\left( \omega _{r}^{2}-{{\omega }^{2}} \right)}{{{\left( \omega _{r}^{2}-{{\omega }^{2}} \right)}^{2}}+{{\left( 2{{\zeta }_{r}}{{\omega }_{r}}\omega \right)}^{2}}}, \\ Y\left( \omega \right)=\sum \limits_{r=1}^{NDOF} \frac{\omega \sum \limits_{s=1}^{NELP} {{{\hat{P}}}_{sr}}\sum \limits_{s=1}^{NELP} {{{\hat{P}}}_{rs}}\left( 2{{\zeta }_{r}}{{\omega }_{r}}\omega \right)}{{{\left( \omega _{r}^{2}-{{\omega }^{2}} \right)}^{2}}+{{\left( 2{{\zeta }_{r}}{{\omega }_{r}}\omega \right)}^{2}}}. \\ \end{array}\end{eqnarray} \tag{ 60 }$$

It should be noted that the optimal load resistance can be substituted back into equation (58) to give optimal power harvesting. Moreover, the multimode FRF representing the transverse displacement relative to the input transverse acceleration can be obtained as

$$\begin{eqnarray}&&\begin{array}{lllllllllllllll} \frac{{{a}_{r}}\left( j\omega \right)}{-{{\omega }^{2}}{{w}_{base}}{{e}^{j\omega \,t}}}\;=-\frac{1}{\omega _{r}^{2}-{{\omega }^{2}}+j2{{\zeta }_{r}}{{\omega }_{r}}\omega } \\ \end{array}\begin{array}{lllllllllllllll} \quad \times (\frac{\left( \sum \limits_{s=1}^{NELP} j\omega {{P}_{Ds}}+{{R}_{L}} \right){{{\hat{Q}}}_{r}}}{\sum \limits_{s=1}^{NELP} j\omega {{P}_{Ds}}-\frac{1}{{{R}_{load}}}-\sum \limits_{r=1}^{NDOF} \frac{j\omega \sum \limits_{s=1}^{NELP} {{{\hat{P}}}_{sr}}\sum \limits_{s=1}^{NELP} {{{\hat{P}}}_{rs}}}{\omega _{r}^{2}-{{\omega }^{2}}+j2{{\zeta }_{r}}{{\omega }_{r}}\omega }}). \\ \end{array}\end{eqnarray} \tag{ 61 }$$

Before equation (61) is modified into the transverse displacement FRF with respect to the input transverse acceleration at any position along the unimorph beam (x), the characteristic transverse motion of the unimorph beam can be reformulated in terms of the natural normal mode multiplied by the time-dependent transverse displacement, using equations (21) and (46) to give

$$\begin{eqnarray}\begin{array}{lllllllllllllll} w\left( x,t \right) & = & \mathop{{{{\bf \Phi }}^{{\rm e}}}\left( x \right)}\limits_{\left( 1\times 4 \right)}\mathop{{{{\boldsymbol{ \varphi }} }^{{\rm e}}}}\limits_{\left( 4\times r \right)}\mathop{{\boldsymbol{a}} \left( t \right)}\limits_{\left( r\times 1 \right)} \\ {} & = & \mathop{{{{\boldsymbol{ \lambda }} }^{e}}\left( x \right)}\limits_{\left( 1\times r \right)}\mathop{{\boldsymbol{a}} \left( t \right)}\limits_{\left( r\times 1 \right)} \\ {} & = & \sum \limits_{{\rm r}=1}^{{\rm NDOF}} {{\lambda }^{e}}_{r}\left( x \right){{a}_{r}}\left( t \right). \\ \end{array}\end{eqnarray} \tag{ 62 }$$

Moreover, the characteristic motions for tip mass offset over the interval $L\leqslant x\leqslant L+{{L}_{tip}}$ can be formulated using the characteristic motions of the tip of the beam. As mentioned previously, note that since the tip mass offset is assumed to undergo rigid body motion, its motion depends solely upon the tip of the beam's transverse rectilinear and rotation or slope motions. To this point, the transverse displacement response of the tip mass offset can be formulated as

$$\begin{eqnarray}&&\begin{array}{lllllllllllllll} {\boldsymbol{w}} \left( L+{{L}_{tip}},t \right)={\boldsymbol{w}} \left( L,t \right)+{{L}_{tip}}\frac{{\rm d}{\boldsymbol{w}} \left( L,t \right)}{{\rm d}x} \\ \quad =\;\left( \mathop{{{{\boldsymbol{ \Phi }} }^{{\rm e}}}\left( L \right)}\limits_{\left( 1\times 4 \right)}+{{L}_{tip}}\mathop{\frac{{\rm d}{{{\boldsymbol{ \Phi }} }^{{\rm e}}}\left( L \right)}{{\rm d}x}}\limits_{\left( 1\times 4 \right)} \right)\mathop{{{{\boldsymbol{ \varphi }} }^{{\rm e}}}}\limits_{\left( 4\times r \right)}\mathop{a\left( t \right)}\limits_{\left( r\times 1 \right)} \\ \quad =\;\mathop{{{{\bf \lambda }}^{{\rm e}}}\left( L+{{L}_{tip}} \right)}\limits_{\left( 1\times r \right)}\mathop{{\boldsymbol{a}} \left( t \right)}\limits_{\left( r\times 1 \right)}=\sum \limits_{r=1}^{NDOF} {{{\rm \lambda }}^{e}}_{r}\left( L+{{L}_{tip}} \right){{a}_{r}}\left( t \right). \\ \end{array}\end{eqnarray} \tag{ 63 }$$

Corresponding to equations (61) and (62), the FRF that relates the multimode transverse displacement to the input base acceleration of the unimorph beam can be formulated as

$$\begin{eqnarray}&&\begin{array}{lllllllllllllll} \frac{w\left( x,j\omega \right)}{-{{\omega }^{2}}{{w}_{base}}{{{\rm e}}^{j\omega \,t}}}=-\sum \limits_{r=1}^{NDOF} \left[ \frac{{{{\rm \lambda }}^{e}}_{r}\left( x \right)}{\omega _{r}^{2}-{{\omega }^{2}}+j2{{{\rm \zeta }}_{r}}{{\omega }_{r}}\omega } \right. \\ \quad \times (\frac{\left( \sum \limits_{s=1}^{NELP} j\omega {{P}_{Ds}}+{{R}_{L}} \right){{{\hat{Q}}}_{r}}}{\sum \limits_{s=1}^{NELP} j\omega {{P}_{Ds}}-\frac{1}{{{R}_{load}}}-\sum \limits_{r=1}^{NDOF} \frac{j\vec{\omega }\sum \limits_{s=1}^{NELP} {{{\hat{P}}}_{sr}}\sum \limits_{s=1}^{NELP} {{{\hat{P}}}_{rs}}}{\omega _{r}^{2}-{{{\rm \omega }}^{2}}+j{\rm 2}{{\zeta }_{r}}{{\omega }_{r}}\omega }}){\mkern 1mu} ]. \\ \end{array}\end{eqnarray} \tag{ 64 }$$

For tip mass offset, the multimode transverse displacement FRFs can also be formulated in terms of equations (61) and (63) to give

$$\begin{eqnarray}&&\begin{array}{lllllllllllllll} \frac{w\left( L+{{L}_{tip}},j\omega \right)}{-{{\omega }^{2}}{{w}_{base}}{{e}^{j\omega \,t}}}=-\sum \limits_{r=1}^{NDOF} [\frac{{{\lambda }^{e}}_{r}\left( L+{{L}_{{\rm tip}}} \right)}{\omega _{r}^{2}-{{\omega }^{2}}+j2{{\zeta }_{r}}{{\omega }_{r}}\omega } \\ (\frac{\left( \sum \limits_{s=1}^{NELP} j\omega {{P}_{Ds}}+{{R}_{L}} \right){{{\hat{Q}}}_{r}}}{\sum \limits_{s=1}^{NELP} j\omega {{P}_{Ds}}-\frac{1}{{{R}_{load}}}-\sum \limits_{r=1}^{NDOF} \frac{j\omega \sum \limits_{s=1}^{NELP} {{{\hat{P}}}_{sr}}\sum \limits_{s=1}^{NELP} {{{\hat{P}}}_{rs}}}{\omega _{r}^{2}-{{\omega }^{2}}+j2{{\zeta }_{r}}{{\omega }_{r}}\omega }})]. \\ \end{array}\end{eqnarray} \tag{ 65 }$$

For experimental study, a laser Doppler vibrometer was used to capture dynamic responses from the small beam and tip mass structures under input base motion, thus giving the absolute motion. For this case, equations (64) and (65) need to be modified in terms of the absolute motion. Since the combination of the base and relative motions defines the absolute motion, the multimode FRFs of absolute transverse displacement and velocity at any position along the unimorph beam and the tip mass offset can be formulated as

$$\begin{eqnarray}\begin{array}{lllllllllllllll} \frac{{{w}_{abs}}\left( x,j\omega \right)}{-{{\omega }^{2}}{{w}_{base}}{{e}^{j\omega \,t}}} & = & \frac{{{w}_{base}}{{e}^{j\omega t}}+w\left( x,j\omega \right)}{-{{\omega }^{2}}{{w}_{base}}{{e}^{j\omega t}}} \\ {} & = & -\frac{1}{{{\omega }^{2}}}+\frac{w\left( x,j\omega \right)}{-{{\omega }^{2}}{{w}_{base}}{{e}^{j\omega t}}} \\ \frac{{{{\dot{w}}}_{abs}}\left( x,j\omega \right)}{-{{\omega }^{2}}{{w}_{base}}{{e}^{j\omega t}}} & = & \frac{\frac{{\rm d}}{{\rm d}t}\left[ {{w}_{base}}{{e}^{j\omega t}}+w\left( x,j\omega \right) \right]}{-{{\omega }^{2}}{{w}_{base}}{{e}^{j\omega t}}} \\ {} & = & \frac{1}{j\omega }+j\omega \frac{w\left( x,j\omega \right)}{-{{\omega }^{2}}{{w}_{base}}{{e}^{j\omega t}}} \\ \frac{{{w}_{abs}}\left( L+{{L}_{tip}},j\omega \right)}{-{{\omega }^{2}}{{w}_{base}}{{e}^{j\omega \,t}}} & = & \frac{{{w}_{base}}{{e}^{j\omega t}}+w\left( L+{{L}_{tip}},j\omega \right)}{-{{\omega }^{2}}{{w}_{base}}{{e}^{j\omega t}}} \\ {} & = & -\frac{1}{{{\omega }^{2}}}+\frac{w\left( L+{{L}_{tip}},j\omega \right)}{-{{\omega }^{2}}{{w}_{base}}{{e}^{j\omega t}}} \\ \frac{{{{\dot{w}}}_{abs}}\left( L+{{L}_{tip}},j\omega \right)}{-{{\omega }^{2}}{{w}_{base}}{{e}^{j\omega \,t}}} & = & \frac{\frac{{\rm d}}{{\rm d}t}\left[ {{w}_{base}}{{e}^{j\omega t}}+w\left( L+{{L}_{tip}},j\omega \right) \right]}{-{{\omega }^{2}}{{w}_{base}}{{e}^{j\omega t}}} \\ {} & = & \frac{1}{j\omega }+j\omega \frac{w\left( L+{{L}_{tip}},j\omega \right)}{-{{\omega }^{2}}{{w}_{base}}{{e}^{j\omega t}}}. \\ \end{array}\end{eqnarray} \tag{ 66 }$$

In this section, the multimode FRFs for the velocity, electrical voltage, current and power outputs using the new electromechanical finite element equations are investigated and validated through experimental studies with variable load resistance (resistive shunt circuit). The input base transverse acceleration onto unimorph structures was chosen to be 1 m s⁻². The piezoelectric unimorph structure was made from PZT PSI-5A4E material (Piezo Systems Inc., Woburn, MA), and the rigid body of the offset tip mass was subsequently attached. The geometrical structures of the piezoelectric unimorph with tip mass offset can be seen in figure 4, and details of the tip mass geometry, including formulas, can be seen in appendix D. Note that, since the tip mass was glued on the top surface of the substructure, the extra substructure length, l_t, was assumed to be a rigid body contributing to the tip mass offset. The complete properties of the unimorph are given in table 1. In figure 4, the unimorph length, L, and width, b, with piezoelectric thickness, h_p , and substructure (brass) thickness, h_s , were set to 60 mm, 6 mm, 0.127 mm and 0.5 mm, respectively. The tip mass configurations were calculated according to the geometry and material properties. The dimensions of tip mass offset with l_t , h_t and b (width) were set to 15 mm, 10 mm and 6 mm, respectively.

Zoom In Zoom Out Reset image size

Table 1.  Properties of the piezoelectric unimorph system.

Material properties	Piezoelectric	Brass
Young's modulus, ${{\bar{c}}_{11}}$ (GPa)	66	105
Density, $\rho$ (kg m−3)	7800	9000
Piezoelectric constant, d31 (pm V−1)	−190	—
Permittivity, $\varepsilon _{33}^{T}$ (F m−1)	1800 ${{\varepsilon }_{{\rm o}}}$	—
Permittivity of free space, ${{\varepsilon }_{{\rm o}}}$ (pF m−1)	8.854	—

As shown in the experimental setup in figure 5, the cantilevered piezoelectric unimorph beam with the tip mass offset was clamped to the base structure under input dynamic excitation. The input acceleration from the base structure was measured using a B&K accelerometer. The B&K exciter connected to the power amplifier was used as the input dynamic excitation. The wave function generator connected to the power amplifier was used to control the input excitation using the sine sweep from the frequencies of 5 Hz to 1000 Hz, with frequency increments of 0.3125 Hz and 40 averages controlled from the B&K Fast Fourier Transform (FFT) Pulse analyzer software, in order to obtain a very smooth frequency response and to avoid aliasing components. Dynamic velocity of the tip mass offset was measured using a Polytec laser vibrometer with a measurement range full scale of 500 mm s⁻¹ and a low pass filter setting of 22 kHz. Before operating the laser, a small reflecting film was attached to the top of the tip mass. Moreover, all signal measurements from the accelerometer, piezoelectric unimorph, and vibrometer were connected to the sophisticated B&K FFT Analyzer 3560B data acquisition device, which has five channels with a maximum frequency span up to 25.6 kHz. All processing signals through the analyzer displayed the measurement results using the FFT pulse software. Note that all sensitivities for scaling factors given from the B&K accelerometer of 9.633 mV m s⁻², and from the laser vibrometer of 125 mm s⁻¹ V (8 V m s⁻¹), were inserted into the FFT pulse in order to obtain the correct results. Since the acceleration is given in per-unit acceleration in m s⁻², all results of the FRFs given from the experiment and the numerical analysis in the next stage are given in per-unit acceleration in m s⁻². Note that the B&K pulse software also provides the facility (a unit organizer) for modifying the units of amplitude of FRF; the amplitudes from the laser measurement can be stated in the velocity or displacement per-unit acceleration.

Zoom In Zoom Out Reset image size

The first four modes of the unimorph beam, as shown in figure 6, can be used to identify the smoothness of mode shapes using the different sizes of a number of discretized elements. Although the resonant frequencies with discretization using a lower number of elements seem to be quite close to each other, the mode shapes seem to be quite coarse. These results show that the proposed 50 electromechanical discretized elements appear to be sufficient to provide accurate results that can give smooth mode shapes, including the natural frequencies. With the 50 electromechanical discretized elements, the strain node along the unimorph beam, as shown in figure 7, can be determined using the inflection point of the second order differential shape function of the strain–displacement relationship, as shown in equation (48). As the surfaces of the piezoelectric material are covered with the electrode layers, the inflection point can be used to identify the location of strain nodes along the beam, indicating the change of strains (tension and compressive strains) between two segments of the piezoelectric structure. This can avoid the reduction of the power output at the second or higher modes. Further details of the parametric strain node cases associated with the experimental study, in terms of electrical connection patterns, can be discussed in future investigations. Our main purpose in this numerical study is to present the new electromechanical finite element vibration techniques and show the validation with experimental studies.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

In figure 8(a), the absolute tip velocity FRFs and experimental results for the first mode show very good agreement under the variable load resistance. As we can see, higher amplitudes can be achieved when load resistances and frequencies shift from short circuits to open circuits. By viewing figure 8(b), the maximum tip absolute velocity responses at the short circuit and open circuit resonance frequencies of 18.5 Hz and 18.9 Hz were reached for load resistances approaching the lower and the higher values (from short circuits to open circuit load resistances), respectively. This behaviour can also be seen more obviously in figure 8(c). In many cases of vibration, the maximum dynamic displacement is generally avoided for most structures. It should be noted here that the power harvesting response can be shown to give the best results without having the maximum displacement, as we will discuss further in the next section. Note that the resonance frequency and amplitude response can shift due to resistive shunt damping from the variable load resistance, where this affects the electromechanical behaviour of the piezoelectric element. In addition, the physical behaviour of the piezoelectric unimorph, involving piezoelectric couplings and internal capacitance, also adds electromechanical damping and electromechanical stiffness. At this point, we realized that the damping effects consist of the mechanical and electrical components, due to the electromechanical behaviour of the dimensional structure, external load resistance and material properties of the piezoelectric structure and substructure. Identification of mechanical damping with very low load resistance can be obtained by matching the amplitude of the experimental and theoretical tip absolute velocity and voltage FRFs. The mechanical damping ratio of the first mode, ${{\zeta }_{1}}=0.0162$, under a very low load resistance of 562 Ω was obtained (approaching the short circuit resistance) because using the actual short circuit load resistance (R_load  = 0) will lead the theoretical voltage FRF to be zero, where the tip absolute displacement cannot be identified. This situation cannot be used to identify the voltage and displacement FRF behaviour of the electromechanical system.

Zoom In Zoom Out Reset image size

With variation in the load resistance, the first mode of the voltage FRFs also varies while increasing load resistance, as shown in figure 9(a). In this case, the resonance frequency shifts to the higher value from short circuit to open circuit load resistances, followed by an increase in the amplitude. Note that when the voltage FRF with load resistance is close to open circuit, the amplitude remains constant.

Zoom In Zoom Out Reset image size

Figure 9(b) shows the trend of voltage behaviour under variable load resistance giving good agreement between numerical and experimental results. In general form, the voltage under short and open circuit resonances indicates a slight increase with increasing load resistance. However, the maximum open circuit resonance amplitude gives a higher value compared to the short circuit resonance when load resistance passes over the transitional or critical point of the short circuit and open circuit amplitudes. Note that the critical point of the amplitude response between short circuit and open circuit resonance frequencies shows the same amplitude. The change of frequency due to the variable load resistance can be seen more clearly in figure 9(c). For electrical current FRFs, as shown in figure 10(a), the trend shows the opposite behaviour of that of the voltage FRFs. The electrical current shows the highest amplitude when the load resistance approaches the open circuit. Again, in the opposite trend to the voltage response, the variation of electric current with variable load resistance under short circuit and open circuit resonance frequencies, as shown in figure 10(b), reaches the maximum amplitude with decreasing load resistance. This trend can be seen more obviously in the 3D current FRF shown in figure 10(c). Another important aspect related to the FRFs is the relationship between dynamic velocity, electrical voltage and current. The absolute velocity amplitude with the load resistance approaching the short circuit seems to increase, whereas the voltage amplitude decreases with increasing electrical current. Conversely, the velocity amplitude with the load resistance approaching open circuit tends to increase, whereas the voltage amplitude increases with decreasing electrical current.

Zoom In Zoom Out Reset image size

However, the power harvester FRFs, especially near the first mode response, as shown in figure 11(a), seem to behave with different trends where the short circuit and open circuit load resistances tend to decrease the amplitude. It can be said that the highest velocity amplitudes at the short circuit and open circuit load resistances do not show the highest power amplitudes. In different configurations with velocity, voltage and current, the variations in power harvesting under short circuit and open circuit resonance frequencies, as seen in figure 11(b), show two different optimum load resistances, giving the optimum amplitudes. In the 3D graph shown in figure 11(c), the optimum power harvesting amplitudes can also be seen clearly at the regions away from the short circuit and open circuit load resistances.

Zoom In Zoom Out Reset image size

Optimal power harvesting, based on the optimal load resistance, can also be seen in figure 12(a), where the absolute maximum power amplitude point coincides with the load resistance of 227 ${\rm k}\Omega$ at the resonance frequency of 18.70 Hz. However, the optimal load resistance in the off-resonance regions seem to approach and overlap with the load resistances of 155 ${\rm k}\Omega$ and 332 ${\rm k}\Omega$, where the experimental results are also very close to the off-resonance regions, except for the driving frequency below 8 Hz due to unwanted noise. By scrutinizing the optimal power harvesting response, the numerical modelling result seen in figure 12(b) shows that ignoring the rotary inertia of the tip mass results in the shift of the resonance frequency to a higher value. The resonance difference is also more pronounced if the model ignores the offset distance of the tip mass measured from its centroid to the end of the unimorph beam length. Neglecting both the offset and rotary inertia of the tip mass results in a larger modelling error. It can be seen that when designing either a micro- or mesoscale power harvesting device, the mathematical model should include the effect of the offset and rotary inertia of the tip mass.

Zoom In Zoom Out Reset image size

The new numerical study proposed in this paper addresses five important points.

• The electromechanical discretization was numerically developed, since the actual piezoelectric material covered by thin electrode layers only generates a single voltage across the electrical load. The electromechanical discretization introduced here consists of the mechanical discretized elements for elemental displacement fields and the parallel electrical discretized elements for multi-output electrical current. The previous published papers related to the piezoelectric finite element have not considered this in their mathematical modelling, although many of them showed different applications.
• Instead of developing the kinematic equations of the unimorph beam, the tip mass offset was formulated, since it contributes significantly to the mass matrices and the input dynamic force vector. The benefit of locating the tip mass with its centroid away from the end of the piezoelectric beam is that it allows for the modelling of various applications, especially MEMS power harvesters. The improved designs avoid contact between the brittle piezoelectric material and the tip mass. The previous published papers related to MEMS have ignored this in their mathematical modelling.
• The numerical techniques provide physical interrelationship between the electrode effect, coupled field equations, kinematic equations, electromechanical discretization and electromechanical dynamic equations.
• The reduced constitutive electromechanical finite element equations were further derived into global normalized scalar form of the electromechanical dynamic equations, in order to formulate the strain mode shape, eigenmode shape and multimode electromechanical FRFs. Unlike analytical techniques, the proposed new numerical solution techniques provided the benefit of analyzing the structure with different scalabilities, including MEMS devices.
• The comparisons between the numerical solution techniques and experimental studies have been discussed in terms of the variable load resistance, giving good agreement.

Overall, the experimental verification and parametric numerical case studies have been explored and have shown close agreement. For this paper, the proposed new electromechanical finite element modelling has proved that the velocity, current and power FRFs have shown very similar trends with the previous published analytical literature. These numerical results were also validated with experimental studies. Moreover, since the trend of optimal power harvesting FRF using different cases of the tip mass was presented, the new method showed that neglecting the rotary inertia of the tip mass resulted in a resonance frequency error. Notably, modelling errors occurred from ignoring the tip mass offset and ignoring both tip mass offset and rotary inertia. The other new aspect shown in using these numerical techniques is that the system models configured normalized strain node and eigenmode shapes using discretized elements for convergence studies. Further details of the frequency analysis using the strain node effect in terms of the electrical circuit connection patterns will be a focus of future studies.

The new numerical techniques can be a very useful tool for analyzing FRFs, displacement mode shape and strain node forms. These techniques can be applied to model different geometrical aspects for laminated structures and MEMS power harvester devices. In this case, since our main concern is presenting the new contribution of numerical techniques, it is important to validate the results obtained from numerical FRFs using experimental studies where the trends of the results were also shown to be very similar with facts established in previous analytical literature. Further applications for new aspects of power harvesting using the proposed numerical techniques can be demonstrated in future research studies.

In figure A1, the location of the asymmetric neutral axis, measured from the y-axis to the top surface of the piezoelectric layer, can be determined using the resultant force balance in the cross section of the unimorph structure to give

$$\begin{eqnarray}&&\begin{array}{l} \mathop \sum \limits_{i=1} {{F}^{(i)}}=\int \limits_{-{{h}_{s}}+\left( {{z}_{n}}-{{h}_{p}} \right)}^{{{z}_{n}}-{{h}_{p}}} \int \limits_{0}^{b} c_{11}^{\left( 1 \right)}z{\rm d}y{\rm d}z \\ \qquad \qquad +\int \limits_{{{z}_{n}}-{{h}_{p}}}^{{{z}_{n}}} \int \limits_{0}^{b} c_{11}^{\left( 2 \right)}z{\rm d}y{\rm d}z \\ \qquad \quad =\,0. \\ \end{array}\end{eqnarray} \tag{ A.1 }$$

Zoom In Zoom Out Reset image size

Manipulating equation (A1), the location of the asymmetric neutral axis can be reduced as

$$\begin{eqnarray}&&{{z}_{n}}=\frac{c_{11}^{\left( 1 \right)}h_{s}^{2}+c_{11}^{\left( 2 \right)}h_{p}^{2}+2c_{11}^{\left( 1 \right)}{{h}_{s}}{{h}_{p}}}{2\left( c_{11}^{\left( 1 \right)}{{h}_{s}}+c_{11}^{\left( 2 \right)}{{h}_{p}} \right)},\end{eqnarray} \tag{ A.2 }$$

where coefficients $\bar{c}_{11}^{\left( 1 \right)}$ and $\bar{c}_{11}^{E}=\bar{c}_{11}^{\left( 2 \right)}$ represent substructure elastic stiffness and piezoelectric elastic stiffness at a constant electric field, respectively.

Corresponding to equations (7) and (11), the total kinetic energy of the system consisting of the elemental unimorph and tip mass offset, as shown in figure 3, can be formulated as

$$\begin{eqnarray}&&KE=K{{E}^{unimorph}}+K{{E}^{tip}},\end{eqnarray} \tag{ B.1 }$$

where

$$\begin{eqnarray}&&\begin{array}{lllllllllllllll} K{{E}^{unimorph}}=\frac{1}{2}\mathop \int \limits_{{{m}^{\left( 1 \right)}}} {{{{\boldsymbol{\dot{R}}} }}^{{\boldsymbol{p}} {{{\boldsymbol{p}} }^{\prime\prime }}}}\left( x,z,t \right)\cdot {{{{\boldsymbol{\dot{R}}} }}^{{\boldsymbol{p}} {{{\boldsymbol{p}} }^{\prime\prime }}}}\left( x,z,t \right){\rm d}{{m}^{\left( 1 \right)}} \\ \quad \quad \quad +\frac{1}{2}\mathop \int \limits_{{{m}^{\left( 2 \right)}}} {{{{\boldsymbol{\dot{R}}} }}^{{\boldsymbol{p}} {{{\boldsymbol{p}} }^{\prime\prime }}}}\left( x,z,t \right)\cdot {{{{\boldsymbol{\dot{R}}} }}^{{\boldsymbol{p}} {{{\boldsymbol{p}} }^{\prime\prime }}}}\left( x,z,t \right){\rm d}{{m}^{\left( 2 \right)}} \\ K{{E}^{tip}}=\frac{1}{2}\mathop \int \limits_{{{m}^{tip}}} {{{{\boldsymbol{\dot{R}}} }}^{m{{m}^{\prime\prime }}}}\left( {{x}_{e+{\rm 1}}},z,t \right)\cdot {{{{\boldsymbol{\dot{R}}} }}^{m{{m}^{\prime\prime }}}}\left( {{x}_{e+{\rm 1}}},z,t \right){\rm d}{{m}^{tip}}. \\ \end{array}\end{eqnarray} \tag{ B.2 }$$

As mentioned previously, the translation and angular velocities of tip mass at the end of the unimorph beam length, L, as shown in figure 1, can be transformed into the local position of the end of the element, ${{x}_{e+1}}$, as shown in figure 3. The first part of equation (B.2) can be expanded as

$$\begin{eqnarray}&&\begin{array}{lllllllllllllll} K{{E}^{unimorph}} \\ =\frac{1}{2}\mathop \int \limits_{{{m}^{\left( 1 \right)}}} \left( {{{{\boldsymbol{\dot{w}}} }}_{base}}\left( t \right)\cdot {{{{\boldsymbol{\dot{w}}} }}_{base}}\left( t \right)+2{{{{\boldsymbol{\dot{w}}} }}_{base}}\left( t \right)\cdot {\boldsymbol{\dot{w}}} \left( x,t \right) \right. \\ \left. +{{z}^{2}}{\boldsymbol{\dot{ \theta }}} \left( x,t \right)\cdot {\boldsymbol{\dot{ \theta }}} \left( x,t \right)+{\boldsymbol{\dot{w}}} \left( x,t \right)\cdot {\boldsymbol{\dot{w}}} \left( x,t \right) \right)\;{\rm d}{{m}^{\left( 1 \right)}} \\ +\frac{1}{2}\mathop \int \limits_{{{m}^{\left( 2 \right)}}} \left( {{{{\boldsymbol{\dot{w}}} }}_{base}}\left( t \right)\cdot {{{{\boldsymbol{\dot{w}}} }}_{base}}\left( t \right)+2{{{{\boldsymbol{\dot{w}}} }}_{base}}\left( t \right)\cdot {\boldsymbol{\dot{w}}} \left( x,t \right) \right. \\ \left. +{{z}^{2}}{\boldsymbol{\dot{ \theta }}} \left( x,t \right)\cdot {\boldsymbol{\dot{ \theta }}} \left( x,t \right)+{\boldsymbol{\dot{w}}} \left( x,t \right)\cdot {\boldsymbol{\dot{w}}} \left( x,t \right) \right)\;{\rm d}{{m}^{\left( 2 \right)}} \\ \end{array}\end{eqnarray} \tag{ B.3 }$$

Note that the first terms for each parenthesis from equation (B.3) must be zero, due to the Lagrangian functional form of the system only depending upon the relative translation and velocity of the unimorph, as expressed in section 4.2. The second terms for each parenthesis contribute to the input dynamic forces, and these terms can be moved into the nonconservative work section, as given in detail in appendix C. For the tip mass kinetic energy, the second part of equation (B.2) can be expanded as

$$\begin{eqnarray}&&\begin{array}{lllllllllllllll} K{{E}^{tip}}=\frac{1}{2}\mathop \int \limits_{{{m}^{tip}}} \left( {{{{\boldsymbol{\dot{w}}} }}_{base}}\left( t \right)\cdot {{{{\boldsymbol{\dot{w}}} }}_{base}}\left( t \right)+2{{x}_{m}}{\boldsymbol{\dot{ \theta }}} \left( {{x}_{e+1}},t \right)\cdot {{{{\boldsymbol{\dot{w}}} }}_{base}}\left( t \right) \right. \\ \ \;+2{{x}_{m}}{\boldsymbol{\dot{ \theta }}} \left( {{x}_{e+1}},t \right)\cdot {\boldsymbol{\dot{w}}} \left( {{x}_{e+1}},t \right)+2{{x}_{c}}{{x}_{m}}{\boldsymbol{\dot{ \theta }}} \left( {{x}_{e+1}},t \right)\cdot {\boldsymbol{\dot{ \theta }}} \left( {{x}_{e+1}},t \right) \\ \ \;+2{{z}_{c}}{{z}_{m}}{\boldsymbol{\dot{ \theta }}} \left( {{x}_{e+1}},t \right)\cdot {\boldsymbol{\dot{ \theta }}} \left( {{x}_{e+1}},t \right)+2{{x}_{c}}{{{{\boldsymbol{\dot{w}}} }}_{base}}\left( t \right)\cdot {\boldsymbol{\dot{ \theta }}} \left( {{x}_{e+1}},t \right) \\ \ \;+2{{{{\boldsymbol{\dot{w}}} }}_{base}}\left( t \right)\cdot {\boldsymbol{\dot{w}}} \left( {{x}_{e+1}},t \right)+{\boldsymbol{\dot{w}}} \left( {{x}_{e+1}},t \right){\boldsymbol{\dot{w}}} \left( {{x}_{e+1}},t \right) \\ \ \;+2{{x}_{c}}{\boldsymbol{\dot{ \theta }}} \left( {{x}_{e+1}},t \right)\cdot {\boldsymbol{\dot{w}}} \left( {{x}_{e+1}},t \right)+x_{c}^{2}{\boldsymbol{\dot{ \theta }}} \left( {{x}_{e+1}},t \right)\cdot {\boldsymbol{\dot{ \theta }}} \left( {{x}_{e+1}},t \right) \\ \ \;+z_{c}^{2}{\boldsymbol{\dot{ \theta }}} \left( {{x}_{e+1}},t \right)\cdot {\boldsymbol{\dot{ \theta }}} \left( {{x}_{e+1}},t \right)+z_{m}^{2}{\boldsymbol{\dot{ \theta }}} \left( {{x}_{e+1}},t \right)\cdot {\boldsymbol{\dot{ \theta }}} \left( {{x}_{e+1}},t \right) \\ \left. \ \;+x_{m}^{2}{\boldsymbol{\dot{ \theta }}} \left( {{x}_{e+1}},t \right)\cdot {\boldsymbol{\dot{ \theta }}} \left( {{x}_{e+1}},t \right) \right)\;{\rm d}{{m}^{tip}}. \\ \end{array}\end{eqnarray} \tag{ B.4 }$$

Again, the first term of equation (B.4) must be zero. Moreover, the second to fifth terms must be zero because they contain the first mass moment of inertias for the tip mass body about the centroid at point g, as shown in figure 3, since the centroid of the tip mass is located on that centroid itself [57]. Again, the sixth and seventh terms can be moved into the nonconservative work section.

Equations (B.3) and (B.4) can simply be rewritten in matrix form as

$$\begin{eqnarray}\begin{array}{lllllllllllllll} K{{E}^{unimorph}} & = & \frac{1}{2}\int \limits_{{{x}_{e}}}^{{{x}_{e+1}}} \mathop \int \limits_{{{A}^{\left( 1 \right)}}} {{z}^{2}}{{\rho }^{\left( 1 \right)}}{\boldsymbol{\dot{ \theta }}} {{\left( x,t \right)}^{T}}{\boldsymbol{\dot{ \theta }}} \left( x,t \right){\rm d}{{A}^{\left( 1 \right)}}{\rm d}x \\ {} & {} & +\frac{1}{2}\int \limits_{{{x}_{e}}}^{{{x}_{e+1}}} \mathop \int \limits_{{{A}^{\left( 2 \right)}}} {{z}^{2}}{{\rho }^{\left( 2 \right)}}{\boldsymbol{\dot{ \theta }}} {{\left( x,t \right)}^{T}}{\boldsymbol{\dot{ \theta }}} \left( x,t \right){\rm d}{{A}^{\left( 2 \right)}}{\rm d}x \\ {} & {} & +\frac{1}{2}\int \limits_{{{x}_{e}}}^{{{x}_{e+1}}} \mathop \int \limits_{{{A}^{\left( 1 \right)}}} {{\rho }^{\left( 1 \right)}}{\boldsymbol{\dot{w}}} {{\left( x,t \right)}^{T}}{\boldsymbol{\dot{w}}} \left( x,t \right){\rm d}{{A}^{\left( 1 \right)}}{\rm d}x \\ {} & {} & +\frac{1}{2}\int \limits_{{{x}_{e}}}^{{{x}_{e+1}}} \mathop \int \limits_{{{A}^{\left( 2 \right)}}} {{\rho }^{\left( 2 \right)}}{\boldsymbol{\dot{w}}} {{\left( x,t \right)}^{T}}{\boldsymbol{\dot{w}}} \left( x,t \right){\rm d}{{A}^{\left( 2 \right)}}{\rm d}x \\ K{{E}^{tip}} & = & I_{0}^{tip}{{x}_{c}}{\boldsymbol{w}} {{\left( {{x}_{e+1}},t \right)}^{T}}{\boldsymbol{\dot{ \theta }}} \left( {{x}_{e+1}},t \right) \\ {} & {} & +\frac{1}{2}I_{0}^{tip}{\boldsymbol{\dot{w}}} {{\left( {{x}_{e+1}},t \right)}^{T}}{\boldsymbol{\dot{w}}} \left( {{x}_{e+1}},t \right) \\ {} & {} & +\frac{1}{2}I_{2}^{tip}{\boldsymbol{\dot{ \theta }}} {{\left( {{x}_{e+1}},t \right)}^{T}}{\boldsymbol{\dot{ \theta }}} \left( {{x}_{e+1}},t \right). \\ \end{array}\end{eqnarray} \tag{ B.5 }$$

The zeroth and second mass moment of inertias of the tip mass offset can be stated as

$$\begin{eqnarray}\begin{array}{lllllllllllllll} I_{0}^{tip} & = & \mathop \int \limits_{{{V}^{tip}}} {{\rho }^{tip}}{\rm d}{{V}^{tip}}, \\ I_{2}^{tip} & = & \mathop \int \limits_{{{V}^{tip}}} {{\rho }^{tip}}\left( z_{m}^{2}+x_{m}^{2} \right)d{{V}^{tip}}+\left( z_{c}^{2}+x_{c}^{2} \right){{m}^{tip}} \\ {} & = & I_{cg}^{tip}+{{m}^{tip}}\left( z_{c}^{2}+x_{c}^{2} \right), \\ \end{array}\end{eqnarray} \tag{ B.6 }$$

where ${\rm d}{{V}^{tip}}={\rm d}{{x}_{m}}{\rm d}{{y}_{m}}{\rm d}{{z}_{m}}$ and definite integral forms of f(y, z) implied from equation (B.5) can be defined as

$$\begin{eqnarray}\begin{array}{lllllllllllllll} \mathop \int \limits_{{{A}^{\left( 1 \right)}}} f\left( y,z \right)\;{\rm d}{{A}^{\left( 1 \right)}} & = & \int \limits_{-{{h}_{s}}+\left( {{z}_{n}}-{{h}_{p}} \right)}^{{{z}_{n}}-{{h}_{p}}} \int \limits_{0}^{b} f\left( y,z \right)\;{\rm d}y{\rm d}z, \\ \mathop \int \limits_{{{A}^{\left( 2 \right)}}} f\left( y,z \right){\rm d}{{A}^{\left( 2 \right)}} & = & \int \limits_{{{z}_{n}}-{{h}_{p}}}^{{{z}_{n}}} \int \limits_{0}^{b} f\left( y,z \right){\rm d}y{\rm d}z, \\ {} & {} & f\left( y,z \right)\in \left\{ 1,z,{{z}^{2}} \right\}. \\ \end{array}\end{eqnarray} \tag{ B.7 }$$

Note that equation (B.7) can be used for any operational integrations, not only in the kinetic energy derivation, but also for other Lagrangian energy forms given from equations (36)–(39).

As mentioned in appendix B, the contributions of nonconservative forces due to the input base excitation on the system were reduced from the kinetic energy. As a result, the virtual work done by nonconservative forces can simply be formulated in terms of their time-dependent response, using the Hamiltonian principle as

$$\begin{eqnarray}&&\delta W{{F}^{base}}=\int \limits_{{{t}_{1}}}^{{{t}_{2}}} {\boldsymbol{F}} .\delta {\boldsymbol{r}} {\rm d}t.\end{eqnarray} \tag{ C.1 }$$

Note that from a mathematical standpoint, the functional form of Lagrangian's principle is reduced from the Hamiltonian principle using the variational method. Therefore, corresponding to equation (39), the virtual nonconservative works can be expanded by using the second terms for each parenthesis from equation (B.3) and the sixth and seventh terms from equation (B.4), as

$$\begin{eqnarray}\begin{array}{lllllllllllllll} \delta W{{F}^{base}} & = & \frac{1}{2}\left( \int \limits_{{{t}_{1}}}^{{{t}_{2}}} \mathop \int \limits_{{{m}^{(1)}}} 2{{{{\boldsymbol{\dot{w}}} }}_{base}}\left( t \right)\cdot \delta {\boldsymbol{\dot{w}}} \left( x,t \right){\rm d}{{m}^{\left( 1 \right)}}{\rm d}t \right. \\ {} & {} & +\int \limits_{{{t}_{1}}}^{{{t}_{2}}} \mathop \int \limits_{{{m}^{(2)}}} 2{{{{\boldsymbol{\dot{w}}} }}_{base}}\left( t \right)\cdot \delta {\boldsymbol{\dot{w}}} \left( x,t \right){\rm d}{{m}^{\left( 2 \right)}}{\rm d}t \\ {} & {} & +\int \limits_{{{t}_{1}}}^{{{t}_{2}}} \mathop \int \limits_{{{m}^{tip}}} 2{{x}_{c}}{{{{\boldsymbol{\dot{w}}} }}_{base}}\left( t \right)\cdot \delta {\boldsymbol{\dot{ \theta }}} \left( {{x}_{e+1}},t \right){\rm d}{{m}^{tip}}{\rm d}t \\ {} & {} & \left. +\int \limits_{{{t}_{1}}}^{{{t}_{2}}} \mathop \int \limits_{{{m}^{tip}}} 2{{{{\boldsymbol{\dot{w}}} }}_{base}}\left( t \right)\cdot \delta {\boldsymbol{\dot{w}}} \left( {{x}_{e+1}},t \right){\rm d}{{m}^{tip}}{\rm d}t \right). \\ \end{array}\end{eqnarray} \tag{ C.2 }$$

Equation (C.2) needs to be further formulated by applying partial integration, and the result can be stated in the matrix form as

$$\begin{eqnarray}\begin{array}{lllllllllllllll} \delta W{{F}^{base}} & = & -\int \limits_{{{t}_{1}}}^{{{t}_{2}}} \int \limits_{{{x}_{e}}}^{{{x}_{e+1}}} \mathop \int \limits_{{{A}^{\left( 1 \right)}}} {{\rho }^{\left( 1 \right)}}{{{{\boldsymbol{\ddot{w}}} }}_{base}}\left( t \right)\delta {\boldsymbol{w}} {{\left( x,t \right)}^{T}}{\rm d}{{A}^{\left( 1 \right)}}{\rm d}x{\rm d}t \\ {} & {} & -\int \limits_{{{t}_{1}}}^{{{t}_{2}}} \int \limits_{{{x}_{e}}}^{{{x}_{e+1}}} \mathop \int \limits_{{{A}^{\left( 2 \right)}}} {{\rho }^{\left( 2 \right)}}{{{{\boldsymbol{\ddot{w}}} }}_{base}}\left( t \right)\delta {\boldsymbol{w}} {{\left( x,t \right)}^{T}}{\rm d}{{A}^{\left( 2 \right)}}{\rm d}x{\rm d}t \\ {} & {} & -\int \limits_{{{t}_{1}}}^{{{t}_{2}}} I_{0}^{tip}{{x}_{c}}{{{{\boldsymbol{\ddot{w}}} }}_{base}}\left( t \right)\delta {\boldsymbol{ \theta }} {{\left( {{x}_{e+1}},t \right)}^{T}}{\rm d}t \\ {} & {} & -\int \limits_{{{t}_{1}}}^{{{t}_{2}}} I_{0}^{tip}{{{{\boldsymbol{\ddot{w}}} }}_{base}}\left( t \right)\delta {\boldsymbol{w}} {{\left( {{x}_{e+1}},t \right)}^{T}}{\rm d}t. \\ \end{array}\end{eqnarray} \tag{ C.3 }$$

Corresponding to the parameters given by equation (B.6), the mass moment of inertias of the tip mass offset, as shown in figure 4, can be formulated. Note that as mentioned previously, the extra unimorph beam length, as shown in figure D1, also contributed to the tip mass offset. In this case, the zeroth mass moment of inertia can be stated as

$$\begin{eqnarray}&&I_{{\rm 0}}^{tip}={{\rho }^{tip}}b{{l}_{t}}{{h}_{t}}+{{\rho }^{\left( 1 \right)}}b{{l}_{t}}{{h}_{s}},\end{eqnarray} \tag{ D.1 }$$

and the second mass moment of inertia of tip mass offset at the point d, located at the end of the unimorph beam, can be reduced to give

$$\begin{eqnarray}\begin{array}{lllllllllllllll} I_{2}^{tip} & = & {{\rho }^{tip}}b{{l}_{t}}{{h}_{t}}\left( \frac{\left( l_{t}^{2}+h_{t}^{2} \right)}{12}+{{\left( \frac{{{h}_{t}}}{2}-{{h}_{p}}-{{{\bar{z}}}_{c}}+{{z}_{n}} \right)}^{2}} \right) \\ {} & {} & +{{\rho }^{\left( 1 \right)}}b{{l}_{t}}{{h}_{s}}\left( \frac{\left( l_{t}^{2}+h_{s}^{2} \right)}{12}+{{\left( \frac{{{h}_{s}}}{2}+{{h}_{p}}-{{z}_{n}}+{{{\bar{z}}}_{c}} \right)}^{2}} \right) \\ {} & {} & +I_{0}^{tip}\left( \bar{z}_{c}^{2}+\bar{x}_{c}^{2} \right), \\ \end{array}\end{eqnarray} \tag{ D.2 }$$

where the offset distances measured from the tip mass centroid to the point d in the x- and z-axes can be formulated, respectively, as

$$\begin{eqnarray}&&{{\bar{x}}_{c}}=\frac{{{\rho }^{{\rm tip}}}bl_{t}^{2}{{h}_{t}}+{{\rho }^{\left( 1 \right)}}bl_{t}^{2}{{h}_{s}}}{2\left( {{\rho }^{tip}}b{{l}_{t}}{{h}_{t}}+{{\rho }^{\left( 1 \right)}}b{{l}_{t}}{{h}_{s}} \right)}\end{eqnarray} \tag{ D.3 }$$

$$\begin{eqnarray}&&\begin{array}{lllllllllllllll} {{{\bar{z}}}_{c}}= \\ \frac{\left[ \begin{array}{lllllllllllllll} {{\rho }^{tip}}b{{l}_{t}}{{h}_{t}}\left( \frac{{{h}_{t}}}{2}+{{z}_{n}}-{{h}_{p}} \right)+{{\rho }^{\left( 1 \right)}}b{{l}_{t}}\left( {{z}_{n}}-{{h}_{p}} \right) \\ \times \left( \frac{{{z}_{n}}-{{h}_{p}}}{2} \right)-{{\rho }^{\left( 1 \right)}}b{{l}_{t}}\left( {{h}_{s}}+{{h}_{p}}-{{z}_{n}} \right)\left( \frac{{{h}_{s}}+{{h}_{p}}-{{z}_{n}}}{2} \right) \\ \end{array} \right]}{{{\rho }^{tip}}b{{l}_{t}}{{h}_{t}}+{{\rho }^{\left( 1 \right)}}b{{l}_{t}}\left( {{z}_{n}}-{{h}_{p}} \right)+{{\rho }^{\left( 1 \right)}}b{{l}_{t}}\left( {{h}_{s}}+{{h}_{p}}-{{z}_{n}} \right)}. \\ \end{array}\end{eqnarray} \tag{ D.4 }$$

Zoom In Zoom Out Reset image size

Alternatively, and with the same result, the second mass moment of inertia of tip mass offset at the point d, located at the end of the unimorph beam, can straightforwardly be formulated as

$$\begin{eqnarray}&&\begin{array}{lllllllllllllll} I_{2}^{tip}={{\rho }^{tip}}b{{l}_{t}}{{h}_{t}}\left( \frac{\left( l_{t}^{2}+h_{t}^{2} \right)}{12}+{{\left( {{z}_{n}}-{{h}_{p}}+\frac{{{h}_{t}}}{2} \right)}^{2}}+{{\left( \frac{{{l}_{t}}}{2} \right)}^{2}} \right) \\ \ \;+{{\rho }^{\left( 1 \right)}}b{{l}_{t}}{{h}_{s}}\left( \frac{\left( l_{t}^{2}+h_{s}^{2} \right)}{12}+{{\left( \frac{{{l}_{t}}}{2} \right)}^{2}}+{{\left( \frac{{{h}_{s}}}{2}-\left( {{z}_{n}}-{{h}_{p}} \right) \right)}^{2}} \right). \\ \end{array}\end{eqnarray} \tag{ D.5 }$$