## 2 Mathematical Model
### Non-Uniform Euler-Bernoulli Beam Equation and Discretization for Undamped Dynamics Analysis  
The proposed unimorph piezoelectric energy harvester includes beam with free geometry design and piezoelectric film on the upper layer of cantilever beam substructure; Figure 1. L is the length of the beam, b is the width of the beam, \(h(x)\) is the thickness of the beam along \(x\) axis. The transverse displacement of the beam at an arbitrary point x and along the neutral axis is denoted by \(v(x,t)\). Using the Newtonian approach, the governing equation of transversal motion for undamped free vibrations of a non-uniform thickness Euler-Bernoulli beam can be expressed by [31],  
\[\rho A(x)\frac{\partial^{2}v(x,t)}{\partial t^{2}}+\frac{\partial^{2}}{ \partial x^{2}}\left[EI(x)\frac{\partial^{2}v(x,t)}{\partial x^{2}}\right]=0, \tag{1}\]  
\(t\) is the time variable, \(\rho\) is the mass density, \(A(x)\) is the cross-sectional area, E is the Young's modulus, and \(I(x)\) is the second moment of area. Boundary conditions of the non-uniform thickness cantilever beam can be defined as,  
\[ \begin{split}[v(x,t)]_{x=0}&=0,\\ \left[\frac{\partial v(x,t)}{\partial x}\right]_{x=0}& =0,\\ \left[EI(x)\frac{\partial^{2}v(x,t)}{\partial x^{2}}\right]_{x= L}&=0,\\ \left[\frac{\partial}{\partial x}\left(EI(x)\frac{\partial^{2}v( x,t)}{\partial x^{2}}\right)\right]_{x=L}&=0\end{split}  \tag{2}\]  
The formal eigenvalue equation for beam is  
\[EI(x)\frac{\partial^{4}y}{\partial x^{4}}=\zeta\rho A(x)y, \tag{3}\]  
or using special function \(\Phi(x)\)  
\[EI(x)\frac{\partial^{4}\Phi(x)}{\partial x^{4}}=\zeta\rho A(x)\Phi(x), \tag{4}\]  
Figure 1: Unimorph cantilevered beam structure: \(L\)—length, \(b\)—width, \(b_{p}\)—thickness of piezoelectric layer, \(h(x)\)—thickness of substructure, \(\bar{V}\)—acceleration of base mechanical excitation, \(R_{\text{load}}\)—electrical load. where \(\zeta\) eigenvalue. For rectangular cross section bean with thickness h(x) and width b, (4) can be written as  
\[E\frac{bh^{3}(x)}{12}\frac{\partial^{4}\Phi(x)}{\partial x^{4}}=\zeta\rho bh(x) \Phi(x), \tag{5}\]  
The 2D simplified Bernoulli-Euler beam element is presented in Figure 2.  
Using Hermite shape functions and nodal point vector for the beam element, the transversal displacement of the beam at a point x is approximated [32]  
\[v(x)= \frac{v_{1}}{L^{3}}\big{(}2x^{3}-3 Lx^{2}+L^{3}\big{)}-\frac{\theta_{1}}{L ^{3}}\big{(}2x^{3}-3 Lx^{2}\big{)}\ \ +\frac{v_{2}}{L^{2}}\Big{(}x^{3}-2 Lx^{2}+L^{2}x\Big{)}+\frac{\theta_{2}}{L^{2}} \Big{(}x^{3}-Lx^{2}\Big{)},\] (6) where \[u^{e}=[v_{1}\ \ \theta_{1}\ v_{2}\ \ \theta_{2}]^{T} \tag{7}\]  
\(u^{e}\) is nodal points displacement vector of the 2D beam element. By calculating the strain and kinetic energy due to the deformation of the beam, [33; 34] the 2D beam-element stiffness K\({}^{e}\) and mass \(M^{e}\) matrices are described. \[K^{e}=\frac{E}{L^{3}}\left[ \begin{array}{cccc}12I&6 LI&-12I&6 LI\\ &4L^{2}I&-6 LI&2L^{2}I\\ \\ \textit{symmetric}&12I&-6 LI\\ \\ &&&4L^{2}I\end{array} \right] \tag{8}\]  
\[M^{e}=\frac{\rho hL}{420}\left[ \begin{array}{cccc}156&22L&54&-13L\\ &4L^{2}&13L&-3L^{2}\\ \\ \textit{symmetric}&156&-22L\\ \\ &&&4L^{2}\end{array} \right] \tag{9}\]  
The finite element discretization leads to the following discretized equation of motion of the undamped free vibration for the global beam system is  
\[[M(h)]\big{\{}\tilde{u}\big{\}}+[K(h)]\big{\{}u\big{\}}=0 \tag{10}\]  
Or after (10), by rearranging it, it takes the following eigenvalue problem form:  
\[[K(h)]\big{\{}\psi\big{\}}-\omega^{2}[M(h)]\big{\{}\psi\big{\}}=0 \tag{11}\]  
Figure 2: Degrees of freedom of the 2D simplified Bernoulli–Euler bending beam element: \(L\)—length, \(v_{1},v_{2}\)—deflection, \(\theta_{1},\theta_{2}\)—slope at the node. where \(\lambda=\omega^{2}\) is representing vibration resonant frequencies (eigenvalues) while \(\psi\) is representing vibration resonant modes (eigenvectors).