# Vibration Analysis of Piezoelectric Cantilever Beams with Bimodular Functionally-Graded Properties
## 4 Numerical Simulation
### Constitutive Equation of Piezoelectrical Materials  
First, we need to describe the input of the piezoelectrical materials parameters in the software. The piezoelectrical materials model in ABAQUS follows the e-form constitutive equation, such that  
\[\left\{ \begin{array}{l}\sigma_{ij}=c_{ijkl}^{E}\varepsilon_{kl}-e_{kij}E_{k }\\ D_{i}=e_{ijkkl}+\lambda_{ik}^{\varepsilon}E_{k}\end{array} \right., \tag{19}\]  
where the stress component is denoted by \(\sigma_{ij}\); the strain component is denoted by \(\varepsilon_{ij}\); the electrical displacement component is denoted by \(D_{i}\); the electrical field strength are denoted by \(E_{k}\); the stiffness coefficient matrix is denoted by \(c_{ijkl}^{E}\); the dielectric constant matrix is denoted by \(\lambda_{ik}^{\varepsilon}\). In piezoelectrical materials, there is a polarization direction which corresponds to \(z\) direction in \(x\)-\(y\)-\(z\) coordinate system (i.e., 3-direction in matrix). In ABAQUS, we need to transform the flexibility coefficient matrix into the stiffness coefficient matrix, such that  
\[\left[s_{ij}\right]=\left[ \begin{array}{cccccc}s_{11}&s_{12}&s_{13}&0&0&0 \\ s_{21}&s_{22}&s_{23}&0&0&0\\ s_{31}&s_{32}&s_{33}&0&0&0\\ 0&0&0&s_{66}&0&0\\ 0&0&0&0&s_{44}&0\\ 0&0&0&0&0&s_{44}\end{array} \right]\Rightarrow\left[c_{ij}\right]=\left[ \begin{array}{cccccc}c_{11}&c_{12}&c_{13}&0&0&0\\ c_{21}&c_{22}&c_{23}&0&0&0\\ c_{31}&c_{32}&c_{33}&0&0&0\\ 0&0&0&c_{66}&0&0\\ 0&0&0&0&c_{44}&0\\ 0&0&0&0&0&c_{44}\end{array} \right]. \tag{20}\]  
Note that the above definition for the variation form of functionally-graded materials and the flexibility coefficient \(s_{ij}=s_{ij}^{0}e^{\alpha_{i}z/h}\) (where \(i=1,2\) due to the bimodular effect), the stiffness coefficient will thus vary with \(c_{ij}=c_{ij}^{0}e^{-\alpha_{i}z/h}\), otherwise \(s_{ij}\) and \(c_{ij}\) cannot satisfy \([s_{ij}][c_{ij}]=[s_{ij}^{0}]e^{\alpha_{i}z/h}[c_{ij}^{0}]e^{-\alpha_{i}z/h}=[E]\) (here \([E]\) is an unit matrix). The piezoelectrical strain constants should, at the same time, be \(e_{ij}=d_{ij}c_{ij}^{E}=d_{ij}^{0}e^{\alpha_{i2}/h}c_{ij}^{0}e^{-\alpha_{i2}/h}=e_{ ij}^{0}\), where \(c^{E}\) is the constant matrix of elastic stiffness in the state of short circuit; thus, the piezoelectrical stress constants matrix is  
\[\left[e_{ij}\right]=\left[e_{ij}^{0}\right]=\left[ \begin{array}{cccccc}0&0&0& 0&e_{15}&0\\ 0&0&0&0&0&e_{15}\\ e_{31}&e_{31}&e_{33}&0&0&0\end{array} \right]. \tag{21}\]  
The constitutive equation for piezoelectric materials is, in the form of matrix, expressed as  
\[\left[ \begin{array}{c}\sigma_{11}\\ \sigma_{22}\\ \sigma_{33}\\ \sigma_{12}\\ \sigma_{13}\\ \sigma_{23}\end{array} \right]=\left[ \begin{array}{cccccc}c_{11}&c_{12}&c_{ 13}&0&0&0\\ c_{12}&c_{11}&c_{13}&0&0&0\\ c_{13}&c_{13}&c_{33}&0&0&0\\ 0&0&0&c_{66}&0&0\\ 0&0&0&0&c_{44}&0\\ 0&0&0&0&0&c_{44}\end{array} \right]\left[ \begin{array}{c}\varepsilon_{11}\\ \varepsilon_{22}\\ \varepsilon_{33}\\ \gamma_{12}\\ \gamma_{13}\\ \gamma_{23}\end{array} \right]-\left[ \begin{array}{cc}0&0&e_{31}\\ 0&0&e_{33}\\ 0&0&e_{33}\\ 0&0&0\\ e_{15}&0&0\end{array} \right]\left[ \begin{array}{c}E_{1}\\ E_{2}\\ E_{3}\end{array} \right] \tag{22}\]  
and  
\[\left[ \begin{array}{c}D_{1}\\ D_{2}\\ D_{3}\end{array} \right]=\left[ \begin{array}{cccccc}0&0&0&0&e_{15}&0\\ 0&0&0&0&e_{15}\\ e_{31}&e_{31}&e_{33}&0&0&0\end{array} \right]\left[ \begin{array}{c}\varepsilon _{11}\\ \varepsilon_{22}\\ \varepsilon_{33}\\ \gamma_{12}\\ \gamma_{13}\\ \gamma_{23}\end{array} \right]+\left[ \begin{array}{cc}\lambda_{11}&0&0\\ 0&\lambda_{11}&0\\ 0&0&\lambda_{33}\end{array} \right]\left[ \begin{array}{c}E_{1}\\ E_{2}\\ E_{3}\end{array} \right] \tag{23}\]  
According to Voigt notation, the vector components (1, 2, 3, 4, 5 and 6) correspond to the second-order tensor mark of double-subscript (11, 22, 33, 13, 23 and 12), respectively. Thus, the above two equations can be written as,  
\[\left[ \begin{array}{c}\sigma_{11}\\ \sigma_{22}\\ \sigma_{33}\\ \sigma_{12}\\ \sigma_{13}\\ \sigma_{23}\end{array} \right]=\left[