Flexoelectric effect on thickness-shear vibration of a rectangular piezoelectric crystal plate

During the past decades, piezoelectric materials have been poured into engineering industries to fabricate various electromechanical transducers based on the physical phenomenon, piezoelectricity, which is capable of mutually converting mechanical and electric energies. The fabricated piezoelectric devices can be classified into actuators and sensors which utilize different functioning modes, such as flexural mode, thickness mode, extension mode, shear mode, twisting mode and radial mode [1]. In communication industry, piezoelectric Quartz crystal resonator based on acoustic technology is an essential frequency control component and provides frequency reference signal with high stability. The functioning mode of crystal resonator is the thickness-shear mode (TSh) which is also a widely used vibration mode for piezoelectric devices, such as surface acoustic wave (SAW) resonators, film bulk (FBAR) resonators, gyroscope and other acoustic wave sensors. Study of TSh vibration of piezoelectric structures is of great importance for their design and operation in engineering applications.

Numerous works have been carried out for the TSh vibration analysis of piezoelectric crystals by Mindlin, Tiersten, Lee and other researchers [2–4] for electroded or unelectroded resonators [5, 6]. To improve the performance of resonators, TSh vibration of piezoelectric plate with nonuniform thickness [7, 8] was investigated and contoured or beveled resonators [4, 9] were developed based on the energy trapping effect [10, 11]. The nonlinear coupling between TSh mode and thickness-stretch mode has also been studied for a rotated Y-cut Quartz resonator considering both material and kinematic nonlinearities [12]. The nonlinear TSh vibration of quartz crystal plates under an electric field was studied by Wu et al [13]. Yang et al [14] studied the middle plane curvature effect on TSh mode of crystal resonator by two-dimensional equations. TSh vibrations of piezoelectric cylindrical shell were also studied by Yang et al [15, 16]. Zhao [17] studied the TSh and thickness-twist modes in an AT-cut Quartz acoustic wave filter. The effect of a liquid layer on TSh vibrations of crystal were investigated by Lee and Jing with closed form solutions [18, 19]. Besides, the acceleration effect on the vibration frequency of TSh mode of an infinite isotropic plate was also investigated by Wu et al [20]. However, in recent years, the advent of the 5G era leads to a large number of electronic components moving towards miniaturization. Engineers made the device smaller and thinner to meet the requirements of higher frequency communication. The thickness of resonant devices become much smaller and the frequencies of the devices become much larger. Therefore, as the main functioning mode of the crystal resonator, TSh vibration is also affected by the smaller thickness of resonator devices. Therefore, it is essential to study the size effect on vibration frequencies of TSh mode in order to accurately design high frequency devices.

Recent studies have found that when the device size is at micro/nano scale, dielectric material exhibits another remarkable electromechanical coupling characteristic, that is, the flexoelectric effect, in which the electric polarization is proportional to the strain gradient [21–23]. Flexoelectric effect is an electromechanical coupling phenomenon widely existing in various dielectric materials and has inspired great research interests by scientists. When the material or structure is at macro size, the polarization caused by the flexoelectric effect is often ignored. However, when the size of the material or structure is at micro/nano size, the influence of the flexoelectric effect becomes very significant, and it becomes larger and larger as the size of the structure decreases [24–26]. Thus, flexoelectric effect has obvious size effect and recent studies have found that flexoelectricity affects the propagation of Lamb wave and Raleigh wave in half-space solids [27–29]. Sanjeev et al [30] investigated the wave velocity change due to boundary types for bedded piezo-structure with flexoelectric effect. Jiao et al [31] investigated the wave propagation through a flexoelectric piezoelectric slab and discussed the results of the influences of two characteristic lengths and the flexoelectric coefficients on the wave propagation. Zhu et al [32] recently investigated the thickness-twist waves in the nanoplates with flexoelectricity and found that flexoelectric effect cannot be neglected to predict wave propagation behavior in nanoscale structures.

As the main functioning mode of piezoelectric resonators, TSh vibration is also supposed to be significantly affected for small size dielectric materials when considering the flexoelectric effect. The thicknesses of commonly used resonators are from dozens of micron to hundreds of micron range and their frequencies are from a few megahertz to thousand megahertz range. Developing an analytical model for TSh vibration considering flexoelectricity and investigating the flexoelectric effect on vibration frequencies of TSh mode are very essential for the design and miniaturization of piezoelectric resonators. At present, there is a lack of theoretical model for the TSh vibration with flexoelectricity, although there are numerous works focusing on the TSh vibration of piezoelectric devices. Thus, in this paper we propose an analytical model for the TSh vibration of piezoelectric crystal considering flexoelectricity and present TSh vibration results of a rectangular piezoelectric crystal plate. The effect of flexoelectricity on the TSh vibration will be detailedly investigated in the following section

Consider a rectangular piezoelectric crystal plate in figure 1 with the length equals 2a and thickness equals 2h. x₁ and x₂ are the longitudinal direction and thickness direction, respectively. We start the modeling with energy density function. Within the assumption of infinitesimal deformation, the general expression for the electric Gibbs free energy density function H can be written as the following equation for a piezoelectric material with flexoelectricity [33].

$$\begin{eqnarray}&&H\left(E,S,\lambda \right)=-\displaystyle \frac{1}{2}{\varepsilon }_{ij}{E}_{i}{E}_{j}+\displaystyle \frac{1}{2}{c}_{ijkl}{S}_{ij}{S}_{kl}-{e}_{ijk}{E}_{k}{S}_{ij}-{\mu }_{ijkl}{E}_{i}{\lambda }_{jkl}\end{eqnarray} \tag{ 1 }$$

where ε , E , c , S , e , μ and λ are the second-order dielectric constant, electric field, fourth-order elastic constant, strain tensor, piezoelectric tensor, fourth-order flexoelectric tensor and strain gradient tensor, respectively.

Zoom In Zoom Out Reset image size

The constitutive equations for the dielectric materials can be expressed in terms of the Gibbs free energy as

$$\begin{eqnarray}&&\begin{array}{l}{T}_{ij}=\displaystyle \frac{\partial H}{\partial {S}_{ij}}={c}_{ijkl}{S}_{kl}-{e}_{ijk}{E}_{k}\\ {\tau }_{jkl}=\displaystyle \frac{\partial H}{\partial {\lambda }_{jkl}}=-{\mu }_{ijkl}{E}_{i}\\ {D}_{i}=-\displaystyle \frac{\partial H}{\partial {E}_{i}}={\varepsilon }_{ij}{E}_{j}+{e}_{ijk}{S}_{jk}+{\mu }_{ijkl}{\lambda }_{jkl}\end{array}\end{eqnarray} \tag{ 2 }$$

where T , τ and D represent the Cauchy stress tensor, higher-order stress and electric displacement vector, respectively.

The elastic strain and strain gradient are defined in terms of displacement u as follows.

$$\begin{eqnarray}&&\begin{array}{l}{S}_{ij}=\displaystyle \frac{1}{2}\left({u}_{i,j}+{u}_{j,i}\right)\\ {\lambda }_{ijk}={S}_{ij,k}=\displaystyle \frac{1}{2}{\left({u}_{i,j}+{u}_{j,i}\right)}_{,k}\end{array}\end{eqnarray} \tag{ 3 }$$

The electric field can be written in terms of electrostatic potential ϕ as follows.

$$\begin{eqnarray}&&{E}_{i}=-{\phi }_{,i}\end{eqnarray} \tag{ 4 }$$

The governing equations and boundary conditions of linear piezoelectricity can be derived from the variational principle. We consider the following variational functional [34]:

$$\begin{eqnarray}\begin{array}{lll}{\rm{\Pi }}\left(u,\phi \right) & = & \displaystyle {\int }_{0}^{{t}_{0}}dt\displaystyle {\int }_{V}\left[\displaystyle \frac{1}{2}\rho {\dot{u}}_{i}{\dot{u}}_{i}+\rho {f}_{i}{u}_{i}-H\left(E,S,\lambda \right)-{\rho }_{e}\phi \right]dV\\ & & +\,\displaystyle {\int }_{0}^{{t}_{0}}dt\displaystyle {\int }_{S}{\bar{t}}_{i}{u}_{i}dS-\displaystyle {\int }_{0}^{{t}_{0}}dt\displaystyle {\int }_{S}{\bar{\sigma }}_{e}\phi dS\end{array}\end{eqnarray} \tag{ 5 }$$

where $\rho ,$ ${f}_{i},$ ${\rho }_{e},$ ${\bar{t}}_{i}$ and ${\bar{\sigma }}_{e}$ are the density, body force, volume charge density, prescribed traction force and free charge density per unit surface area.

The first variation of the above functional Π is

$$\begin{eqnarray}\begin{array}{lll}\delta {\rm{\Pi }} & = & \displaystyle {\int }_{0}^{{t}_{0}}dt\displaystyle {\int }_{V}\left\{\left[{\left({T}_{ij}-{\tau }_{ijk,k}\right)}_{,j}+\rho {f}_{i}-\rho {\ddot{u}}_{i}\right]\delta {u}_{i}+\left({D}_{i,i}-{\rho }_{e}\right)\delta \phi \right\}dV\\ & & -\,\displaystyle {\int }_{0}^{{t}_{0}}dt\displaystyle {\int }_{S}\left[\left({T}_{ij}-{\tau }_{ijk,k}\right){n}_{j}-{\bar{t}}_{i}\right]\delta {u}_{i}dS-\displaystyle {\int }_{0}^{{t}_{0}}dt\displaystyle {\int }_{S}{\tau }_{ijk}{n}_{k}\delta {u}_{i,j}dS\\ & & -\,\displaystyle {\int }_{0}^{{t}_{0}}dt\displaystyle {\int }_{S}\left({D}_{i}{n}_{i}+{\bar{\sigma }}_{e}\right)\delta \phi dS=0\end{array}\end{eqnarray} \tag{ 6 }$$

The further derivation of ${\tau }_{ijk}{n}_{k}\delta {u}_{i,j}$ can be found in the reference [35]. In the present work, we adopt the simplified flexoelectric theory and neglect the higher order stress boundaries. Therefore, the stationary condition of Π is

$$\begin{eqnarray*}&&{\left({T}_{ij}-{\tau }_{ijk,k}\right)}_{,j}+\rho {f}_{i}=\rho {\ddot{u}}_{i}\,{\rm{in}}\,V\end{eqnarray*}$$

$$\begin{eqnarray*}&&{D}_{i,i}={\rho }_{e}\,{\rm{in}}\,V\end{eqnarray*}$$

$$\begin{eqnarray}&&\left({T}_{ij}-{\tau }_{ijk,k}\right){n}_{j}={\bar{t}}_{i}\,{\rm{on}}\,S\end{eqnarray} \tag{ 7 }$$

$$\begin{eqnarray*}&&{D}_{i}{n}_{i}=-{\bar{\sigma }}_{e}\,{\rm{on}}\,S\end{eqnarray*}$$

Ignoring the body force and volume charge density, the governing equations are reduced to

$$\begin{eqnarray}&&\begin{array}{l}{\left({T}_{ij}-{\tau }_{ijk,k}\right)}_{,j}=\rho {\ddot{u}}_{i},\left(\left|{x}_{1}\right|\lt a,\left|{x}_{2}\right|\lt h\right)\\ {D}_{i,i}=0,\left(\left|{x}_{1}\right|\lt a,\left|{x}_{2}\right|\lt h\right),\end{array}\end{eqnarray} \tag{ 8 }$$

For the boundary conditions, if there is no traction force and surface charge, they can be written as the following equations for a finite rectangular plate.

$$\begin{eqnarray}&&\begin{array}{l}{T}_{i2}-{\tau }_{i2k,k}=0,\,{D}_{2}=0,\,\left({x}_{2}=\pm h\right)\\ {T}_{1j}-{\tau }_{1jk,k}=0,\,{D}_{1}=0,\,\left({x}_{1}=\pm a\right)\end{array}\end{eqnarray} \tag{ 9 }$$

For the analysis of TSh vibration, we are mainly interested in the anti-symmetric mode through thickness and symmetric mode in the length so that the displacement field with amplitude A is described by the following equations [34].

$$\begin{eqnarray}&&\begin{array}{l}{u}_{1}=A\,\cos \left(\displaystyle \frac{m\pi {x}_{1}}{2a}\right)\,\sin \,\left(\displaystyle \frac{n\pi {x}_{2}}{2h}\right){e}^{i\omega t},\left(m=1,3,5,\mathrm{...};\,n=1,3,5,\mathrm{...}\right)\\ {u}_{2}=0,\,{u}_{3}=0\end{array}\end{eqnarray} \tag{ 10 }$$

The electric potential function is a two-dimensional function and can be assumed as

$$\begin{eqnarray}&&\phi \left({x}_{1},{x}_{2},t\right)=\phi \left({x}_{1},{x}_{2}\right){e}^{i\omega t}\end{eqnarray} \tag{ 11 }$$

Let $f\left({x}_{1},t\right)=A\,\cos \left(\tfrac{m\pi {x}_{1}}{2a}\right){e}^{i\omega t},$ the nonzero strain components can be calculated as below.

$$\begin{eqnarray}&&{S}_{11}={u}_{1,1}=f^{\prime} \left({x}_{1},t\right)\,\sin \,\left(\displaystyle \frac{n\pi {x}_{2}}{2h}\right),\,\,\,{S}_{12}=\displaystyle \frac{n\pi }{2h}f\left({x}_{1},t\right)\,\cos \,\left(\displaystyle \frac{n\pi {x}_{2}}{2h}\right)\end{eqnarray} \tag{ 12 }$$

And nonzero strain gradients are given below.

$$\begin{eqnarray}&&{S}_{11,1}={u}_{1,11}=f^{\prime\prime} \left({x}_{1},t\right)\,\sin \,\left(\displaystyle \frac{n\pi {x}_{2}}{2h}\right),{S}_{11,2}={u}_{1,12}=\displaystyle \frac{n\pi }{2h}f^{\prime} \left({x}_{1},t\right)\,\cos \,\left(\displaystyle \frac{n\pi {x}_{2}}{2h}\right)\end{eqnarray} \tag{ 13 }$$

$$\begin{eqnarray*}&&{S}_{12,1}={u}_{1,21}=\displaystyle \frac{n\pi }{2h}f^{\prime} \left({x}_{1},t\right)\,\cos \,\left(\displaystyle \frac{n\pi {x}_{2}}{2h}\right),\,\,\,{S}_{12,2}={u}_{1,22}=-\displaystyle \frac{{n}^{2}{\pi }^{2}}{4{h}^{2}}f\left({x}_{1},t\right)\,\sin \,\left(\displaystyle \frac{n\pi {x}_{2}}{2h}\right)\end{eqnarray*}$$

For the present work, we use the properties of piezoelectric Quartz crystal which is 32 point group crystal for modeling and its flexoelectric coefficient matrix can be found in the reference [36] as below.

$$\begin{eqnarray}{\mu }_{6\times 9}=\left[\begin{array}{ccccccccc}{\mu }_{11} & 0 & 0 & 0 & {\mu }_{15} & -{\mu }_{26} & 0 & -{\mu }_{28} & {\mu }_{19}\\ {\mu }_{15} & 0 & 0 & 0 & {\mu }_{11} & {\mu }_{26} & 0 & {\mu }_{28} & {\mu }_{19}\\ {\mu }_{31} & 0 & 0 & 0 & {\mu }_{31} & 0 & 0 & 0 & {\mu }_{39}\\ {\mu }_{41} & 0 & 0 & 0 & -{\mu }_{41} & {\mu }_{46} & 0 & {\mu }_{48} & 0\\ 0 & {\mu }_{52} & {\mu }_{46} & -{\mu }_{52} & 0 & 0 & {\mu }_{48} & 0 & 0\\ 0 & {\mu }_{11}-{\mu }_{15} & -2{\mu }_{46} & {\mu }_{11}-{\mu }_{15} & 0 & 0 & -2{\mu }_{28} & 0 & 0\end{array}\right]\end{eqnarray} \tag{ 14 }$$

The matrix μ reflects the relationship between the induced polarization and applied strain gradient which indicates that it is more suitable to be written as the form of 3 by 18 rather than 6 by 9. With subscript transformation, the flexoelectric coefficient matrix μ becomes the following matrix.

$$\begin{eqnarray}\begin{array}{l}{\mu }_{3\,\times \,18}=\left[\begin{array}{cccccccccc}{\mu }_{1111} & 0 & 0 & {\mu }_{1221} & 0 & 0 & {\mu }_{1331} & 0 & 0 & \\ 0 & {\mu }_{2112} & {\mu }_{2113} & 0 & {\mu }_{2222} & {\mu }_{2223} & 0 & {\mu }_{2332} & 0 & \cdots \\ 0 & {\mu }_{3112} & {\mu }_{3113} & 0 & {\mu }_{3222} & {\mu }_{3223} & 0 & 0 & {\mu }_{3333} & \end{array}\right.\\ \left.\begin{array}{cccccccccc} & {\mu }_{1231} & 0 & 0 & 0 & {\mu }_{1132} & {\mu }_{1133} & 0 & {\mu }_{1122} & {\mu }_{1123}\\ \cdots & 0 & {\mu }_{2232} & {\mu }_{2233} & {\mu }_{2131} & 0 & 0 & {\mu }_{2121} & 0 & 0\\ & 0 & {\mu }_{3232} & 0 & {\mu }_{3131} & 0 & 0 & {\mu }_{3121} & 0 & 0\end{array}\right]\end{array}\end{eqnarray} \tag{ 15 }$$

where

$$\begin{eqnarray*}&&\begin{array}{l}{\mu }_{1111}={\mu }_{2222}={\mu }_{11},\,{\mu }_{2211}={\mu }_{1122}={\mu }_{15},\\ {\mu }_{3311}={\mu }_{3322}={\mu }_{31},\,{\mu }_{2311}={\mu }_{41},\,{\mu }_{2322}=-{\mu }_{41}\\ {\mu }_{1312}={\mu }_{52},\,{\mu }_{1321}=-{\mu }_{52},\,{\mu }_{2323}={\mu }_{1313}={\mu }_{46},\\ {\mu }_{1123}=-{\mu }_{26},\,{\mu }_{2223}={\mu }_{26},\,{\mu }_{1132}=-{\mu }_{28},\,{\mu }_{2232}={\mu }_{28},\\ {\mu }_{1133}={\mu }_{2233}={\mu }_{19},\,{\mu }_{1331}={\mu }_{2332}={\mu }_{48},\,{\mu }_{3333}={\mu }_{39},\\ {\mu }_{1212}={\mu }_{1221}={\mu }_{11}-{\mu }_{15},\,{\mu }_{1213}=-2{\mu }_{46},\,{\mu }_{1231}=-2{\mu }_{48}\end{array}\end{eqnarray*}$$

For arbitrarily cutted Quartz crystal, the flexoelectric coefficient matrix μ has to be transformed by the following relation.

$$\begin{eqnarray}&&{\mu }^{{\prime} }=A\mu {N}^{T}\end{eqnarray} \tag{ 16 }$$

where θ is the cutting angle of crystal, $m=\,\cos \left(\theta \right),\,n=\,\sin \left(\theta \right),$ and

$$\begin{eqnarray*}A=\left[\begin{array}{ccc}1 & 0 & 0\\ 0 & m & n\\ 0 & -n & m\end{array}\right]\end{eqnarray*}$$

$$\begin{eqnarray*}{N}^{T}=\left[\begin{array}{cccccccccccccccccc}1 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0\\ 0 & m & -n & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0\\ 0 & n & m & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0\\ 0 & 0 & 0 & {m}^{2} & 0 & 0 & {n}^{2} & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & -mn & 0 & 0\\ 0 & 0 & 0 & 0 & {m}^{3} & -{m}^{2}n & 0 & m{n}^{2} & -{n}^{3} & 0 & 0 & 0 & 0 & 0 & 0 & 0 & -{m}^{2}n & m{n}^{2}\\ 0 & 0 & 0 & 0 & {m}^{2}n & {m}^{3} & 0 & {n}^{3} & m{n}^{2} & 0 & 0 & 0 & 0 & 0 & 0 & 0 & -m{n}^{2} & -{m}^{2}n\\ 0 & 0 & 0 & {n}^{2} & 0 & 0 & {m}^{2} & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & mn & 0\\ 0 & 0 & 0 & 0 & m{n}^{2} & -{n}^{3} & 0 & {m}^{3} & -{m}^{2}n & 0 & 0 & 0 & 0 & 0 & 0 & 0 & {m}^{2}n & -m{n}^{2}\\ 0 & 0 & 0 & 0 & {n}^{3} & m{n}^{2} & 0 & {m}^{2}n & {m}^{3} & 0 & 0 & 0 & 0 & 0 & 0 & 0 & m{n}^{2} & {m}^{2}n\\ 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & m & 0 & 0 & -n & 0 & 0 & 0 & 0 & 0\\ 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & {m}^{2} & -mn & 0 & -mn & {n}^{2} & 0 & 0 & 0\\ 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & mn & {m}^{2} & 0 & -{n}^{2} & -mn & 0 & 0 & 0\\ 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & n & 0 & 0 & m & 0 & 0 & 0 & 0 & 0\\ 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & mn & -{n}^{2} & 0 & {m}^{2} & -mn & 0 & 0 & 0\\ 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & {n}^{2} & mn & 0 & mn & {m}^{2} & 0 & 0 & 0\\ 0 & 0 & 0 & 2mn & 0 & 0 & -2mn & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & {m}^{2}-{n}^{2} & 0 & 0\\ 0 & 0 & 0 & 0 & 2{m}^{2}n & -2m{n}^{2} & 0 & -2{m}^{2}n & 2m{n}^{2} & 0 & 0 & 0 & 0 & 0 & 0 & 0 & {m}^{3}-m{n}^{2} & {n}^{3}-{m}^{2}n\\ 0 & 0 & 0 & 0 & 2m{n}^{2} & 2{m}^{2}n & 0 & -2m{n}^{2} & -2{m}^{2}n & 0 & 0 & 0 & 0 & 0 & 0 & 0 & {m}^{2}n-{n}^{3} & {m}^{3}-m{n}^{2}\end{array}\right]\end{eqnarray*}$$

The above equations, equations (14)–(16), are introduced in detail in the [36]. For arbitrary θ, the transformed coefficient matrix μ is full of nonzero values. If all strain gradients are taken into consideration, the mathematical model will be very complicate and analytical solution cannot be found by the proposed method. Thus, proper simplification has to be made first. Since the length and width are much larger than the thickness, we can ignore all the strain gradients in x₁ and x₃ directions resulting in two nonzero strain gradients based on the assumed displacement field which are S_11,2 and S_12,2. Secondly, since we focus on the thickness-shear vibration, we are interested in the gradient of shear strain S₁₂ in x₂ direction which plays the most important role of TSh vibration. We simply ignore the gradient of normal strain S₁₁ in x₂ direction. Thus, only the shear strain gradient S_12,2 is left in the present model. Then, two transformed flexoelectric coefficients used in this model can be written as

$$\begin{eqnarray}&&\begin{array}{l}{\mu ^{\prime} }_{1122}=mn{\mu }_{48}+\left({m}^{3}-m{n}^{2}\right){\mu }_{15}-\left({m}^{2}n-{n}^{3}\right){\mu }_{26}\\ {\mu ^{\prime} }_{2122}=-{m}^{3}n{\mu }_{11}+{m}^{2}{n}^{2}{\mu }_{41}-{m}^{2}{n}^{2}{\mu }_{26}-m{n}^{3}{\mu }_{46}\\ +{m}^{3}n{\mu }_{48}+{m}^{2}{n}^{2}{\mu }_{39}\end{array}\end{eqnarray} \tag{ 17 }$$

For simplification, we remove the superscript and use μ₁₁₂₂ and μ₂₁₂₂ instead of ${\mu ^{\prime} }_{1122}$ and ${\mu ^{\prime} }_{2122}$ for deriving the following formulations. From the third equation of constitutive equations, we can get the electric displacements by substituting the nonzero strains and strain gradients as below.

$$\begin{eqnarray}&&\begin{array}{l}{D}_{1}={e}_{11}{u}_{1,1}-{\varepsilon }_{11}{\phi }_{,1}+{\mu }_{1122}{u}_{1,22}\\ {D}_{2}={e}_{26}{u}_{1,2}-{\varepsilon }_{22}{\phi }_{,2}+{\mu }_{2122}{u}_{1,22}\end{array}\end{eqnarray} \tag{ 18 }$$

Substituting the above equations into the second equation of the governing equations, we get

$$\begin{eqnarray}&&\begin{array}{l}{D}_{1,1}={e}_{11}{u}_{1,11}-{\varepsilon }_{11}{\phi }_{,11}+{\mu }_{1122}{u}_{1,122}=0\\ {D}_{2,2}={e}_{26}{u}_{1,22}-{\varepsilon }_{22}{\phi }_{,22}+{\mu }_{2122}{u}_{1,222}=0\end{array}\end{eqnarray} \tag{ 19 }$$

Then, we can solve the following second order derivatives of electric potential function which are expressed in terms of displacement function.

$$\begin{eqnarray}&&\begin{array}{l}{\phi }_{,11}=\displaystyle \frac{{e}_{11}{u}_{1,11}+{\mu }_{1122}{u}_{1,122}}{{\varepsilon }_{11}}\\ {\phi }_{,22}=\displaystyle \frac{{e}_{26}{u}_{1,22}+{\mu }_{2122}{u}_{1,222}}{{\varepsilon }_{22}}\end{array}\end{eqnarray} \tag{ 20 }$$

The integration of equation (20) is performed with respect to x₁ and x₂ as below.

$$\begin{eqnarray}&&\begin{array}{l}{\phi }_{,1}=\displaystyle \frac{{e}_{11}{u}_{1,1}+{\mu }_{1122}{u}_{1,22}}{{\varepsilon }_{11}}+{g}_{1}\left({x}_{2}\right)\\ {\phi }_{,2}=\displaystyle \frac{{e}_{26}{u}_{1,2}+{\mu }_{2122}{u}_{1,22}}{{\varepsilon }_{22}}+{g}_{2}\left({x}_{1}\right)\end{array}\end{eqnarray} \tag{ 21 }$$

In the following, we consider two electric boundary conditions: open circuit condition and short circuit condition.

• (a)  
  Open circuit conditionThe following traction free boundary condition and electric displacement condition have to be satisfied for the open circuit condition [37].

  $$\begin{eqnarray}&&\begin{array}{l}{T}_{12}-{\tau }_{122,2}=0,\,{D}_{1}=0,\,\left({x}_{1}=\pm a\right)\\ {T}_{12}-{\tau }_{122,2}=0,\,{D}_{2}=0,\,\left({x}_{2}=\pm h\right)\end{array}\end{eqnarray} \tag{ 22 }$$

  According to the constitutive equation, we know

  $$\begin{eqnarray}&&\begin{array}{l}{T}_{12}={c}_{66}{u}_{1,2}+{e}_{26}{\phi }_{,2}\\ {\tau }_{122}={\mu }_{1122}{\phi }_{,1}+{\mu }_{2122}{\phi }_{,2}\end{array}\end{eqnarray} \tag{ 23 }$$

  We substitute equations (21) and (23) into (22), leading to the following equation.

  $$\begin{eqnarray}&&\begin{array}{l}\left({c}_{66}+\displaystyle \frac{{e}_{26}^{2}}{{\varepsilon }_{22}}\right){u}_{1,2}+{e}_{26}{g}_{2}\left({x}_{1}\right)-{\mu }_{1122}\left(\displaystyle \frac{{e}_{11}{u}_{1,12}+{\mu }_{1122}{u}_{1,222}}{{\varepsilon }_{11}}+{{g}_{1}}^{^{\prime} }\left({x}_{2}\right)\right)\\ -{\mu }_{2122}\left(\displaystyle \frac{{\mu }_{2122}{u}_{1,222}}{{\varepsilon }_{22}}\right)=0\left({x}_{1}=\pm a\right)\,{\rm{or}}\,\left({x}_{2}=\pm h\right)\end{array}\end{eqnarray} \tag{ 24 }$$

  Substituting the displacement function into the traction free boundary condition at the up and bottom surfaces, we can obtain two equations. Since we assume the anti-symmetric mode through thickness and n is odd value in the displacement field u₁, $\cos \left(\tfrac{n\pi }{2}\right)=0$ and the corresponding terms vanish resulting in the following two equations.

  $$\begin{eqnarray}&&\begin{array}{l}{e}_{26}{g}_{2}\left({x}_{1}\right)-{\mu }_{1122}{{g}_{1}}^{{\prime} }\left(h\right)=0\\ {e}_{26}{g}_{2}\left({x}_{1}\right)-{\mu }_{1122}{{g}_{1}}^{{\prime} }\left(-h\right)=0\end{array}\end{eqnarray} \tag{ 25 }$$

  Then, we can get the following relations.

  $$\begin{eqnarray}&&\begin{array}{l}{{g}_{1}}^{{\prime} }\left(h\right)=g^{\prime} \left(-h\right)\\ {g}_{2}\left({x}_{1}\right)=\displaystyle \frac{{\mu }_{1122}}{{e}_{26}}{{g}_{1}}^{{\prime} }\left(h\right)={B}_{1}=const\end{array}\end{eqnarray} \tag{ 26 }$$

  The electric potential is calculated by integrating the second equation of equation (21) with respect to x₂ and the following expression can be obtained.

  $$\begin{eqnarray}&&\phi =\displaystyle \frac{{e}_{26}{u}_{1}+{\mu }_{2122}{u}_{1,2}}{{\varepsilon }_{22}}+{B}_{1}{x}_{2}+h\left({x}_{1}\right)\end{eqnarray} \tag{ 27 }$$

  where h(x₁) has to be determined by electric displacement condition.Substituting the potential function equations (27) into (18), we get the following electric displacement condition.

  $$\begin{eqnarray}&&\begin{array}{l}{D}_{1}={e}_{11}{u}_{1,1}-{\varepsilon }_{11}\left(\displaystyle \frac{{e}_{26}{u}_{1,1}+{\mu }_{2122}{u}_{1,12}}{{\varepsilon }_{22}}+{h}^{^{\prime} }\left({x}_{1}\right)\right)\\ +{\mu }_{1122}{u}_{1,22}=0,\left({x}_{1}=\pm a\right)\\ {D}_{2}={B}_{1}=0,\,\left({x}_{2}=\pm h\right)\end{array}\end{eqnarray} \tag{ 28 }$$

  Thus, substituting the displacement field u₁ into D₁ at ${x}_{1}=\pm a$ and using the condition that m is odd value, we get

  $$\begin{eqnarray}&&\begin{array}{l}-\displaystyle \frac{m\pi }{2a}{e}_{11}\,\sin \,\left(\displaystyle \frac{m\pi }{2}\right)\,\sin \,\left(\displaystyle \frac{n\pi {x}_{2}}{2h}\right)-\displaystyle \frac{{\varepsilon }_{11}}{{\varepsilon }_{22}}\left(-\displaystyle \frac{m\pi }{2a}{e}_{26}\,\sin \,\left(\displaystyle \frac{m\pi }{2}\right)\,\sin \,\left(\displaystyle \frac{n\pi {x}_{2}}{2h}\right)\right.\\ \left.-\displaystyle \frac{mn{\pi }^{2}}{4bh}{\mu }_{2122}\,\sin \,\left(\displaystyle \frac{m\pi }{2}\right)\cos \left(\displaystyle \frac{n\pi {x}_{2}}{2h}\right)+h\left({x}_{1}\right)\right)=0\\ \displaystyle \frac{m\pi }{2a}{e}_{11}\,\sin \,\left(\displaystyle \frac{m\pi }{2}\right)\,\sin \,\left(\displaystyle \frac{n\pi {x}_{2}}{2h}\right)-\displaystyle \frac{{\varepsilon }_{11}}{{\varepsilon }_{22}}\left(\displaystyle \frac{m\pi }{2a}{e}_{26}\,\sin \,\left(\displaystyle \frac{m\pi }{2}\right)\,\sin \,\left(\displaystyle \frac{n\pi {x}_{2}}{2h}\right)\right.\\ \left.+\,\displaystyle \frac{mn{\pi }^{2}}{4bh}{\mu }_{2122}\,\sin \,\left(\displaystyle \frac{m\pi }{2}\right)\,\cos \,\left(\displaystyle \frac{n\pi {x}_{2}}{2h}\right)+h\left({x}_{1}\right)\right)=0\end{array}\end{eqnarray} \tag{ 29 }$$

  We add the above two equations and easily find that ${h}^{{\prime} }\left({x}_{1}\right)=0.$ Thus, the electric potential for open circuit condition is

  $$\begin{eqnarray}&&\phi =\displaystyle \frac{{e}_{26}{u}_{1}+{\mu }_{2122}{u}_{1,2}}{{\varepsilon }_{22}}\end{eqnarray} \tag{ 30 }$$

  Substituting the displacement function u₁ into the above equation, we obtain the following explicit form of electric potential function.

  $$\begin{eqnarray}&&\begin{array}{l}\phi \left({x}_{1},{x}_{2},t\right)=A\,\cos \left(\displaystyle \frac{m\pi {x}_{1}}{2a}\right)\left(\displaystyle \frac{{e}_{26}}{{\varepsilon }_{22}}\,\sin \left(\displaystyle \frac{n\pi {x}_{2}}{2h}\right)+\displaystyle \frac{n\pi }{2h}\displaystyle \frac{{\mu }_{2122}}{{\varepsilon }_{22}}\,\cos \left(\displaystyle \frac{n\pi {x}_{2}}{2h}\right)\right)\\ {e}^{i\omega t}\left(m=1,3,5,7\cdots ;n=1,3,5,7\cdots \right)\end{array}\end{eqnarray} \tag{ 31 }$$

  Then, we can calculate the frequency from the governing equation by substituting the constitutive equation, potential function and displacement function into equation (8). We simplify it and get the following equation of motion.

  $$\begin{eqnarray}&&{c}_{11}{u}_{1,11}+\displaystyle \frac{{e}_{11}^{2}{u}_{1,11}}{{\varepsilon }_{11}}+{c}_{66}{u}_{1,22}+\displaystyle \frac{{e}_{26}^{2}{u}_{1,22}}{{\varepsilon }_{22}}-\displaystyle \frac{{\mu }_{1122}^{2}{u}_{1,2222}}{{\varepsilon }_{11}}-\displaystyle \frac{{\mu }_{2122}^{2}{u}_{1,2222}}{{\varepsilon }_{22}}=\rho {\ddot{u}}_{1}\end{eqnarray} \tag{ 32 }$$

  Substituting the displacement into equation (32) and simplifying it, we get

  $$\begin{eqnarray}&&\begin{array}{l}{\left(\displaystyle \frac{m\pi }{2a}\right)}^{2}\left({c}_{11}+\displaystyle \frac{{e}_{11}^{2}}{{\varepsilon }_{11}}\right)+{\left(\displaystyle \frac{n\pi }{2h}\right)}^{2}\left({c}_{66}+\displaystyle \frac{{e}_{26}^{2}}{{\varepsilon }_{22}}\right)\\ +{\left(\displaystyle \frac{n\pi }{2h}\right)}^{4}\left(\displaystyle \frac{{\mu }_{1122}^{2}}{{\varepsilon }_{11}}+\displaystyle \frac{{\mu }_{2122}^{2}}{{\varepsilon }_{22}}\right)={\omega }^{2}\rho \end{array}\end{eqnarray} \tag{ 33 }$$

  Finally, the frequency can be obtained and expressed as below.

  $$\begin{eqnarray}&&\begin{array}{l}\omega =\sqrt{\displaystyle \frac{\left({c}_{11}+\tfrac{{e}_{11}^{2}}{{\varepsilon }_{11}}\right){\left(\tfrac{m\pi }{2a}\right)}^{2}+\left({c}_{66}+\tfrac{{e}_{26}^{2}}{{\varepsilon }_{22}}\right){\left(\tfrac{n\pi }{2h}\right)}^{2}+\left(\tfrac{{\mu }_{1122}^{2}}{{\varepsilon }_{11}}+\tfrac{{\mu }_{2122}^{2}}{{\varepsilon }_{22}}\right){\left(\tfrac{n\pi }{2h}\right)}^{4}}{\rho }}\\ \left(m=1,3,5,7\cdots ;n=1,3,5,7\cdots \right)\end{array}\end{eqnarray} \tag{ 34 }$$
• (b)  
  Short circuit condition

For short circuit condition, the traction free boundary condition is same with the open circuit condition, while the electric boundary conditions are not. The boundary conditions are given as below for this case.

$$\begin{eqnarray}&&\begin{array}{l}{T}_{12}-{\tau }_{122,2}=0,\,{D}_{1}=0,\,\left({x}_{1}=\pm a\right)\\ {T}_{12}-{\tau }_{122,2}=0,\,\phi =U{e}^{i\omega t},\,\left({x}_{2}=\pm h\right)\end{array}\end{eqnarray} \tag{ 35 }$$

Same with the open circuit condition, we can get the same potential function equation (27) by applying the traction free boundary condition. Substituting ${x}_{1}=\pm a$ into D₁, we can obtain the result ${h}^{{\prime} }\left({x}_{1}\right)=0$ as same as the previous case. At two surfaces ${x}_{2}=\pm h,$ let U = 0 and the potential boundary condition leads to the following two equations.

$$\begin{eqnarray}&&\begin{array}{l}A\left(\displaystyle \frac{{e}_{26}}{{\varepsilon }_{22}}\,\cos \,\left(\displaystyle \frac{m\pi {x}_{1}}{2a}\right)\,\sin \,\left(\displaystyle \frac{n\pi }{2}\right)+\displaystyle \frac{{\mu }_{2122}}{{\varepsilon }_{22}}\displaystyle \frac{n\pi }{2h}\,\cos \,\left(\displaystyle \frac{m\pi {x}_{1}}{2a}\right)\,\cos \,\left(\displaystyle \frac{n\pi }{2}\right)\right)+{B}_{1}h=0\\ A\left(\displaystyle \frac{{e}_{26}}{{\varepsilon }_{22}}\,\cos \,\left(\displaystyle \frac{m\pi {x}_{1}}{2a}\right)\,\sin \,\left(-\displaystyle \frac{n\pi }{2}\right)+\displaystyle \frac{{\mu }_{2122}}{{\varepsilon }_{22}}\displaystyle \frac{n\pi }{2h}\,\cos \,\left(\displaystyle \frac{m\pi {x}_{1}}{2a}\right)\,\cos \,\left(-\displaystyle \frac{n\pi }{2}\right)\right)-{B}_{1}h=0\end{array}\end{eqnarray} \tag{ 36 }$$

Thus, for odd value of n, we can get

$$\begin{eqnarray}&&{B}_{1}=A{(-1)}^{\displaystyle \frac{n+1}{2}}\displaystyle \frac{{e}_{26}}{{\varepsilon }_{22}h}\,\cos \left(\displaystyle \frac{m\pi {x}_{1}}{2a}\right)\end{eqnarray} \tag{ 37 }$$

Therefore, we can get the potential function for the short circuit condition as below.

$$\begin{eqnarray}&&\begin{array}{l}\phi \left({x}_{1},{x}_{2},t\right)=A\,\cos \left(\displaystyle \frac{m\pi {x}_{1}}{2a}\right)\left(\displaystyle \frac{{e}_{26}}{{\varepsilon }_{22}}\,\sin \left(\displaystyle \frac{n\pi {x}_{2}}{2h}\right)+\displaystyle \frac{n\pi }{2h}\displaystyle \frac{{\mu }_{2122}}{{\varepsilon }_{22}}\,\cos \left(\displaystyle \frac{n\pi {x}_{2}}{2h}\right)\right.\\ \left.+{\left(-1\right)}^{\displaystyle \frac{n+1}{2}}\displaystyle \frac{{e}_{26}}{{\varepsilon }_{22}h}{x}_{2}\right){e}^{i\omega t}\left(m=1,3,5,7\cdots ;n=1,3,5,7\cdots \right)\end{array}\end{eqnarray} \tag{ 38 }$$

The procedure of solving the frequency function is same with the open circuit condition and the frequency function is also same with equation (34). For the sake of brevity, these derivations are not repeated.

As a numerical example, we consier a typical pieozelectric crystal resonator for TSh vibration analysis. For the material properties, we choose the properties of AT-cut Quartz crystal with θ equals 35.25° for numerical analysis [38], where c₁₁ = 86.74 GPa, c₆₆ = 29.01 GPa, e₁₁ = 0.171 C m⁻², e₂₆ = −0.095 C m⁻², ε₁₁ = 39.91 × 10⁻¹² C V⁻¹ m⁻¹, ε₂₂ = 39.82 × 10⁻¹² C V⁻¹m⁻¹, ρ = 2649 kg m⁻³. The reported flexoelectric coefficients for piezoelectric materials such as PZT, BaTiO₃, PVDF and other piezoelectric and dielectric materials varies from hundreds of nC/m to thousands of nC/m [39–41]. However, the flexoelectric coefficients of Quartz crystal materials have not been measured and reported yet. Therefore, based on the flexoelectric coefficients calculated by lattice dynamics and first-principle theory, and experiment measurements of other dielectric materials, we assume that flexoelectric coefficients μ₁₁₂₂ and μ₂₁₂₂ vary from 100 nC m⁻¹ to 300 nC m⁻¹ in this work for investigation of their effect on TSh vibration frequencies.

Figure 2 shows the fundamental frequency ω₁₁ for length to thickness ratio a/h = 5, 10 and 20. It is found that the fundamental frequency decreases dramatically with the increase of thickness value. The length to thickness ratio has slight effect on the fundamental frequency but not very significantly. The thickness range in this graph is commonly used value for resonators with working frequency from several megahertz to one hundred megahertz. In the following, we investigate the flexoelectric effect and the results are shown in figures 3–5 for a/h = 5, 10 and 20, respectively. It should be noted that μ₁₁₂₂ and μ₂₁₂₂ are chosen same in this study. The nondimensional result Ω₁₁ is normalized by ${\omega }_{0}=\sqrt{{c}_{66}/\rho }/4h$ which is the fundamental frequency of TSh mode of an infinite crystal plate without piezoelectricity. From these graphs, it is found that the frequencies are slightly higher than ω₀ with both piezoelectricity and flexoelectricity. The nondimensional frequency is significantly influenced when considering the flexoelectricity, especially for the thickness less than 50 μm. With the increase of flexoelectric coefficients, the nondimensional frequency also increase. However, with the increase of thickness, it is found that the nondimensional frequency decreases and converges to a certain value implying that the flexoelectric effect can be neglected for resonators with large thickness. The length to thickness ratio is also found to be essential for the fundamental frequency, where the nondimensional frequency decreases when a/h ratio increases. Tables 1–3 give the normalized frequencies for n, m = 1, 3, 5, 7, 9 with and without flexoelectric effect for different length to thickness ratios. The results are obtained with μ₁₁₂₂ = μ₂₁₂₂ = 100 nC m⁻¹ and h = 50 μm. From the tables, it should be reminded that the value of nondimensional frequency is slightly changed due to flexoelectricity. For example, it can be seen that the fifth decimal place changes with flexoelectricity for (1, 1) mode, while the fourth and third decimal places change for (3, 1) mode and (5, 1) mode. For other modes, we can get the similar conclusion that frequency change is more significant for higher frequencies due to flexoelectricity. The relative change of the first mode can reach several ppm to dozens of ppm (parts per million). However, this value is remarkable for precise high frequency devices since it is usually specified in terms of ppm and the results are meaningful for design of accurate resonators.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Figure 6 shows the first four TSh modes (u₁) plotted in (x₁, x₂) plane for a/h = 10 and amplitude A = 1. Due to the assumption of anti-symmetry of TSh mode, these modes are anti-symmetric through thickness, but symmetric in x₁ direction. The obtained electric potential functions are also plotted in figures 7 and 8 for open circuit condition and short circuit condition, respectively. Based on the expressions we obtained for electric potential, these two distributions are totally different to meet different electric displacement and potential conditions at the surfaces. The potential contributed by flexoelectricity is also shown in figure 9 which is same for both electric boundary conditions. It is found that the magnitude of electric potential contributed by flexoelectricity is almost 10% of the total electric potential and is much smaller than piezoelectricity for each mode. These potential figures are sufficient to help understand the potential distributions of TSh modes contributed by piezoelectricity and flexoelectricity.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

In this work, we developed a theoretical model for TSh vibration of piezoelectric crystal with flexoelectric effect. The model is derived based on the assumption of anti-symmetric mode through thickness and symmetric mode in length. Exact frequencies and modes of TSh vibration are obtained for a piezoelectric Quartz plate with assumed flexoelectric coefficients. The TSh vibration frequencies are very sensitive to the flexoelectric effect for all modes and significantly affected when the piezoelectric plate is thin. However, the flexoelectric effect can be neglected when the plate is moderate thick in this work. We analyzed the frequency change due to the factors of flexoelectric coefficients, thickness and length to thickness ratio. The results could help understand the flexoelectric effect on the TSh vibration of high frequency piezoelectric resonators with small thickness and help for design of such small size resonators.

This work was supported by the National Natural Science Foundation of China (Grant no. 11702150), Zhejiang Natural Science Foundation (Grant no. LY21A020003), Natural Science Foundation of Ningbo (Grant no. 202003N4015), the project of Key Laboratory of Impact and Safety Engineering (Ningbo University), Ministry of Education (Grant no. CJ202009), the Technology Innovation 2025 Program of Municipality of Ningbo (Grant no. 2019B10122).

The data that support the findings of this study are available upon reasonable request from the authors.