## 2 Mathematical Model
### Piezoelectric Beam Discretization  
The piezoelectric effect enables conversion between mechanical and electrical energy. When a crystal is deformed, it generates electric polarization proportional to the strain, resulting in a potential difference. The equations of piezoelectricity combine the momentum PDE (partial differential equation) with the charge conservation PDE of Electrostatics. PDEs for the mechanical field and for the electrostatic field are coupled by the piezoelectric effect. The mechanical field is described by [PERSON]'s equation  
\[f_{V}+\
abla\cdot[\sigma]=\rho a \tag{26}\]  
where \(f_{V}\) stands for external volume forces and \(\rho\)a is the vector of inertial forces (material density \(\rho\) and acceleration \(a\)). The dot product of the nabla operator \(\
abla\) and the Cauchy stress tensor \([\sigma]\) denotes the divergence of the material stresses. Assuming electrically insulating material and neglecting electromagnetic wave phenomena, [PERSON]'s equations for the electric field reduce to  
\[\
abla\cdot D=q_{\rm e} \tag{27}\]  
\[E=-\
abla V_{\rm e} \tag{28}\]  
with the electric _flux_ density \(D\), the electric charge density \(q_{e}\) on the electrodes, the electric potential \(V_{e}\), and the electric field \(E\), respectively. Equations (26)-(28) are coupled by the constitutive equations of linear piezoelectricity. For the PZT beam it can be formulated with 3-1 mode of piezoelectric constant operation and 3-3 effect of piezoelectric permittivity as,  
\[T_{1}=c_{11}S_{1}-e_{31}^{\top}E_{3} \tag{29}\]  
\[D_{3}=e_{31}S_{1}+e_{33}^{\ast}E_{3} \tag{30}\]  
\(T_{1}\) and \(S_{1}\) represent the vectors of mechanical stress and strain in Voigt notation \(c_{11}\), \(e_{31}\) and \(e_{33}^{\ast}\) are the material tensors for mechanical stiffness at constant electric field, piezoelectric moduli, and electric permeability at constant mechanical strain, respectively. Combining Equations (26)-(30), we arrive at the coupled PDEs for piezoelectric materials in Voigt notation. Coupled field piezoelectric equations in primary dependent variables for mechanical displacement \(u\) and electric potential \(V_{e}\) in piezoelectric medium are as follows:  
\[-\varrho\ddot{u}+\
abla\cdot[c:(\
abla u)]+\
abla\cdot[e\cdot(\
abla v)]+f=0 \tag{31}\]  
\[\
abla\cdot[e:(\
abla u)]-\
abla\cdot[e\cdot(\
abla v)]-q=0 \tag{32}\]  
Or in index notation:  
\[-\varrho\ddot{u}_{i}+c_{11}u+e_{31}\lor+f=0 \tag{33}\]  
\[e_{31}u-e_{33}\lor-q=0 \tag{34}\]  
In the FEM (using Comsol Multiphysics), Equations (33) and (34) are converted into their variational form and the finite element discretized equations for the piezoelectric material are obtained using the standard finite element formulation for the electrostatic equilibrium from Maxwell's equation, the mechanical equilibrium, and the constitutive relation of piezoelectricity together with the elementary and natural boundary conditions as [37]  
\[M_{uu}\ddot{u}+C_{uu}\dot{u}+K_{uu}u+K_{u\varphi}\varphi=R \tag{35}\]  
\[\
abla\cdot D=q_{\rm e} \tag{27}\]0  
where \(M_{uu}\), \(C_{uu}\), \(K_{uu}\), \(K_{u\varphi}\), and \(K_{\varphi\varphi}\) represent, respectively, the mass matrix, mechanical damping matrix, mechanical stiffness matrix, piezoelectric coupling matrix, and dielectric stiffness matrix for the piezoelectric material, respectively, while \(\ddot{u}\), \(\dot{u}\), \(u\), \(R\), \(\varphi\), and \(q\) represent, respectively, the mechanical acceleration vector, the mechanical velocity vector, the mechanical displacement vector, the mechanical external force vector, the electrical potential vector, and the external electric charge vector. Superscript \(T\) stands for the transpose of the matrix. FEM equations of elements can be derived from (35) and (36) and represent dynamic behavior of piezoelectric material in the following form:  
\[ \begin{bmatrix}M_{uu}^{e}&0\\ 0&0\end{bmatrix}\begin{bmatrix}\ddot{u}^{e}\\ \ddot{\varphi}^{e}\end{bmatrix} + \begin{bmatrix}C_{uu}^{e}&0\\ 0&0\end{bmatrix}\begin{bmatrix}\dot{u}^{e}\\ \dot{\varphi}^{e}\end{bmatrix}