## 3 Construction of 3D Piezo Finite Element Model  
(1) Piezoelectric equation  
The coupling of the mechanical and electrical properties of linear piezoelectric materials could be characterized by the following formula (Equations (7) and (8)), and the heat effect could be ignored. \[T=c^{E}S-eE \tag{7}\]  
\[D=eS+\varepsilon^{S}E \tag{8}\]  
where \(T\) is stress, \(c^{E}\) is the elastic matrix under constant strain, \(S\) is the strain, \(e\) is piezoelectric constant, \(E\) is electric field intensity, \(D\) is electric displacement, \(\varepsilon^{S}\) is dielectric constant under the constant electric field. (2) Material property  
Material PZT-4 of the IEEE standard was selected [27], and the parameters were as shown in Equations (9)-(11). \[e=\left[ \begin{array}{cccc}0&0&-5.2\\ 0&0&-5.2\\ 0&0&15.1\\ 0&12.7&0\\ 12.7&0&0\\ 0&0&0\end{array} \right] \tag{9}\]  
\[c^{E}=\left[ \begin{array}{cccccc}1.390&0.778&0.743&0&0&0\\ 0.778&1.390&0.743&0&0&0\\ 0.743&0.743&1.154&0&0&0\\ 0&0&0&0.256&0&0\\ 0&0&0&0&0.256&0\\ 0&0&0&0&0&0.256\end{array} \right] \tag{10}\]  
\[\varepsilon^{S}=\left[ \begin{array}{cccc}6.45&0&0\\ 0&6.45&0\\ 0&0&5.62\end{array} \right] \tag{11}\]  
(3) Element Stiffness matrix  
The hexahedron element, shown in Figure 1, was chosen for the element type of the structure. The displacement field could be expressed by the interpolation function \(N\) and its coordinates, as indicated in Equation (12). \[ \begin{array}{l}x=x_{1}N_{1}+x_{2}N_{2}+x_{3}N_{3}+x_{4}N_{4}\\ y=y_{1}N_{1}+y_{2}N_{2}+y_{3}N_{3}+y_{4}N_{4}\\ z=z_{1}N_{1}+z_{2}N_{2}+z_{3}N_{3}+z_{4}N_{4}\end{array}  \tag{12}\]The interpolation function (or so-called shape function) \(N\) was given by Equation (13) [23]. \[N=\frac{1}{8}\left[ \begin{array}{c}(1-x_{1})(1-x_{2})(1-x_{3})\\ (1+x_{1})(1-x_{2})(1-x_{3})\\ (1+x_{1})(1+x_{2})(1-x_{3})\\ (1-x_{1})(1+x_{2})(1-x_{3})\\ (1-x_{1})(1-x_{2})(1+x_{3})\\ (1+x_{1})(1+x_{2})(1+x_{3})\\ (1-x_{1})(1+x_{2})(1+x_{3})\end{array} \right] \tag{13}\]  
The deduction of the equations of displacement and potential in terms of FEM was discussed. Here, the relationship between \(E\) and \(\phi\) is shown in Equation (14) [28]. \[E=-\overset{\rightarrow}{grad}\phi=-[ \begin{array}{cccc}\overset{ \rightarrow}{grad}N_{1}^{e}&\overset{\rightarrow}{grad}N_{2}^{e}&\ldots& \overset{\rightarrow}{grad}N_{m}^{e}\end{array} ] \tag{14}\]  
This could be rewritten as Equations (15)-(17). \[E=-[B_{\phi}^{e}]\Phi^{e} \tag{15}\]  
\[[B_{\phi}^{e}]=[ \begin{array}{cccc}B_{\phi 1}^{e}&B_{\phi 2}^{e}&\ldots& \ B_{\phi m}^{e}\end{array} ] \tag{16}\]  
\[[B_{\phi i}^{e}]=\left[ \begin{array}{c}\frac{\partial N_{i}^{e}}{\partial N_ {i}^{e}}\\ \frac{\partial N_{i}^{e}}{\partial y}\\ \frac{\partial N_{i}^{e}}{\partial z}\end{array} \right] \tag{17}\]  
where \(\Phi^{e}\) is the potential of the point and \(B_{\phi i}^{e}\) is the derivative of the vector of shape function \(N\). Similarly, as for the displacement field, we used Equations (18)-(20). \[S=[B^{e}]U^{e} \tag{18}\]  
\[[B^{e}]=[ \begin{array}{cccc}B_{1}^{e}&B_{2}^{e}&\ldots&B_{m}^{e}\end{array} ] \tag{19}\]  
Figure 1: 3D hexahedron element and its coordinate system. \[[B_{i}^{e}]=\left[ \begin{array}{ccc}\frac{\partial N_{i}^{e}}{\partial x}&0&0\\ 0&\frac{\partial N_{i}^{e}}{\partial y}&0\\ 0&0&\frac{\partial N_{i}^{e}}{\partial z}\\ 0&\frac{\partial N_{i}^{e}}{\partial z}&\frac{\partial N_{i}^{e}}{\partial y} \\ \frac{\partial N_{i}^{e}}{\partial y}&\frac{\partial N_{i}^{e}}{\partial z}&0 \end{array} \right] \tag{20}\]  
Applying the [PERSON]'s equation [29; 30], the matrix form of the equilibrium equation was obtained as Equation (21). \[\left[ \begin{array}{cc}k_{uu}&k_{u\phi}\\ k_{\phi u}&-k_{\phi\phi}\\ \end{array} \right]\left[ \begin{array}{c}u\\ \phi\\ \end{array} \right]=\left[ \begin{array}{c}F\\ Q\\ \end{array} \right] \tag{21}\]  
The coefficients of the element stiffness matrix were as Equation (22).