# Thickness-Stretch Vibration of an Infinite Piezoelectric Plate with Flexoelectricity
## 2 Mathematical Modeling  
### Three-Dimensional Governing Equations  
For the present work, we start with the general three-dimensional equations. For the piezoelectric plate with the consideration of flexoelectricity, we adopt the Gibbs energy density function as the energy density function for modeling. The Gibbs energy density function \(H\) can be expressed as in Equation (1), which considers the coupling between the strain gradient and electric field [23].  
\[H(E,S,\lambda)=-\frac{1}{2}\epsilon_{ij}E_{i}E_{j}+\frac{1}{2}c_{ijkl}S_{il}S_{ kl}-e_{ijk}E_{k}S_{ij}-\mu_{ijkl}E_{i}\lambda_{jkl} \tag{1}\]  
where \(S\) is the strain tensor, \(\lambda\) is the strain gradient tensor, \(E\) is the electric field, \(c\) is the elastic constant, \(e\) is the piezoelectric constant, \(\varepsilon\) is the dielectric constant, and \(\mu\) is the flexoelectric coefficient.  
The constitutive equations can be derived from the energy density function, as given in Equation (2).  
\[T_{ij}=\frac{\partial H}{\partial S_{ij}}=c_{ijkl}S_{kl}-e_{ijk}E_{k} \tag{2a}\] \[\tau_{jkl}=\frac{\partial H}{\partial\lambda_{jkl}}=-\mu_{ijkl}E_{i} \tag{2b}\]\[D_{i}=-\frac{\partial H}{\partial E_{i}}=\epsilon_{ij}E_{j}+e_{ijk}S_{jk}+\mu_{ ijkl}\lambda_{jkl} \tag{2c}\]  
where \(\mathbf{T}\) is the stress tensor, \(\mathbf{\tau}\) is the higher order stress tensor, and \(\mathbf{D}\) is the electric displacement.  
The linear strains (\(S_{ij}\)) and strain gradients (\(\lambda_{ijk}\)) are calculated by the following geometric relations.  
\[S_{ij}=\frac{1}{2}\left(u_{i,j}+u_{j,j}\right) \tag{3a}\] \[\lambda_{ijk}=S_{ij,k}=\frac{1}{2}\left(u_{i,j}+u_{j,i}\right)_{k} \tag{3b}\]  
The electric field is calculated by the partial differentiation of electric potential, as given in Equation (4).  
\[E_{i}=-\varphi_{j} \tag{4}\]  
The governing equations and boundary conditions of extended linear piezoelectricity can be derived from the variational principle [24], considering the following variational function:  
\[\Pi(u,\varphi)=\int_{0}^{t_{0}}dt\int\limits_{V}\left[\frac{1}{2}\rho u_{i}u_{ i}+\rho f_{i}u_{i}-H(E,S,\lambda)-\rho_{e}\phi\right]dV+\int_{0}^{t_{0}}dt\int \limits_{S}\overline{t}_{i}u_{i}dS-\int_{0}^{t_{0}}dt\int\limits_{S} \overline{\sigma}_{e}\varphi dS \tag{5}\]  
where \(\rho,f_{i}\), \(\rho_{e}\), \(\overline{t}_{i}\) and \(\overline{\sigma}_{e}\) are the density, body force, volume charge density, prescribed traction force, and free charge density per unit surface area.  
Calculating the first variation and the stationary condition gives the following:  
\[\left(T_{ij}-\tau_{ijk,k}\right)_{,j}+\rho f_{i}=\rho\ddot{u}_{i}\text{ in V} \tag{6a}\] \[D_{i,i}=\rho_{e}\qquad\text{ in V}\] (6b) \[\left(T_{ij}-\tau_{ijk,k}\right)n_{j}=\overline{t}_{i}\text{ on S}\] (6c) \[D_{i}n_{i}=-\overline{\sigma}_{e}\text{ on S}\] (6d) \[\tau_{ijk}n_{j}n_{k}=0\text{ on S} \tag{6e}\]  
The governing equations are reduced to the following equations by ignoring the body force and volume charge density:  
\[\left(T_{ij}-\tau_{ijk,k}\right)_{,j}=\rho\ddot{u}_{i} \tag{7a}\] \[D_{i,j}=0 \tag{7b}\]