A Semi-Analytical Solution for the Thickness-Vibration of Centrally Partially-Electroded Circular AT-Cut Quartz Resonators

Vibration frequencies and modes for the thickness-shear vibrations of infinite partially-electroded circular AT-cut quartz plates are obtained by solving the two-dimensional (2D) scalar differential equation derived by Tiersten and Smythe. The Mathieu and modified Mathieu equations are derived from the governing equation using the coordinate transform and the collocation method is used to deal with the boundary conditions. Solutions of the resonant frequencies and trapped modes are validated by those results obtained from COMSOL software. The current study provides a theoretical way for figuring out the vibration analysis of circular quartz resonators.


Introduction
Nowadays, acoustic wave resonators are widely used for frequency generation and control in telecommunications, as well as mass and acceleration sensors. Quartz is the most widely used crystal for resonators, due to its piezoelectric effects. The desire for better resonators has increased with the rapid development of electric devices, which calls for thorough studies on the vibration of the quartz plate. In most applications, quartz crystal plate works with thickness-shear (TS) vibrations [1]. However, because of the complicated properties of quartz crystal, such as anisotropy and electro-mechanical coupling, it is difficult to obtain an exact three-dimensional analytical solution. In order to solve the problem, Mindlin and his co-authors firstly developed a two-dimensional plate theory for solving the vibrations of elastic and piezoelectric plates [2][3][4][5], which is especially suitable for analyzing low-order modes of thickness-shear vibrations. Many researchers have adopted Mindlin's theory to solve vibration of quartz plates in various situations [6][7][8][9][10][11]. Lee [12] established another effective two-dimensional plate theory using trigonometric expansion to express the electric potential and displacement, by which he tried to analyze the stretch and forced vibrations of crystal plates [13,14]. Furthermore, Tiersten and Smythe suggested a two-dimensional scalar differential equation for high-order modes in AT-cut [15] and Stress Compensated-cut (SC-cut) [16] quartz plates, which proved to be quite accurate yet simple. In many of the two-dimensional studies, the plates were assumed to be infinitely large in the x 1 and x 3 directions, and the effects induced by boundaries were considered by in-plane wavenumbers in the x 1 and x 3 directions, which are small compared with the wavenumber in the x 2 direction. This simplification has been shown to be effective in analyzing actual plate vibrations [17][18][19]. In this paper, we will adopt this two-dimensional scalar differential equation for the analysis of various TS modes.
For most quartz resonators, the shape of the crystal plates or the electrodes on the top and bottom of the plates is either rectangular or circular. It is true that rectangular geometry is widely used for Sensors 2017, 17, 1820 2 of 13 general miniaturized quartz resonators. For some specific applications, such as base stations, which do not have strong limitation on space, a circular design of the quartz plate and/or the electrodes on it can show better performance for energy trapping resonance characteristics than a rectangular plate. It is well known that electrode corners in rectangular design can cause electric field concentration, which is associated with degrading effects or even failure of the materials involved. That problem disappears or is very much reduced when a circular design is used, since in this way the electric field concentration can be avoided with a minimal difference in manufacturing. In the case of rectangular plates, the two-dimensional scalar differential equations can be easily decoupled and solved in the Cartesian coordinates, and most of theoretical and numerical studies on TS vibrations mentioned above focused on rectangular plates. In the case of circular resonators, however, the variables in the two-dimensional scalar differential equations cannot be separated easily due to material anisotropy. Therefore, it is necessary to develop a semi-analytical method for circular quartz resonators.
Due to the difficulty of theoretical study on circular resonators, some researchers have turned to numerical solutions. Yong et al. [20] firstly applied Finite Element Method (FEM) to solve the thickness-shear motions of circular crystal plates, followed by Wang et al. [21]. Furthermore, Liu et al. [22,23] developed a new differential quadrature hierarchical FEM method, which can improve accuracy and calculation efficiency. Generally speaking, the FEM approach requires a large amount of time for computation, especially the post process of extracting resonance frequencies and modes is extremely complicated. For theoretical deduction of TS vibration modes, Tiersten proposed a perturbation method, where an incremental item due to anisotropy is added to a nearby isotropic solution [15]. Wang et al. [24] applied the method to solve TS vibration problems of circular plates in order to obtain a truncated 2D equation, while Yang et al. [25] discussed the vibration modes of a circular crystal plate with transversely varying thickness. More recently, He et al. solved the two-dimensional scalar differential equation for a fully-unelectroded finite circular AT-cut quartz plate by coordinate transform and variable separation in elliptical coordinates, by which the original equations were rewritten into the Mathieu and modified Mathieu equations [26]; their article showed better results than the previous perturbation method.
For the actual application of resonators, the electrode usually covers part of the surface on either the top or the bottom, so a partially-electroded model is more accurate in describing an actual resonator than the fully-electroded model. In this paper, the TS vibrations of infinite circular AT-cut quartz plates partially covered with circular electrodes are examined in detail. In Section 2, Tiersten's two-dimensional scalar differential equations are introduced and transformed into elliptical coordinates for decoupling. By introducing boundary and constitutive conditions, the infinite plate's TS vibration frequencies and mode shapes are calculated theoretically, with the numerical models illustrated and discussed in Section 3. In Section 4, FEM simulation by the commercial software COMSOL is described. The FEM solutions are compared with the theoretical results and show good agreement. A thorough understanding of TS vibrations for partially-electroded circular quartz plates can be a theoretical guide for the design, calibration and manufacturing optimization of circular quartz resonators.

Governing Equations
An infinite AT-cut circular quartz plate with two circular identical electrodes placed concentrically on both the top and bottom was studied, as shown in Figure 1. The plate has a diameter of 2R 1 and thickness of 2h; each electrode has a diameter of 2R 0 and thickness of h e . We set the plate radius R 1 = +∞ for simplicity. The vibrations of the electroded region and the unelectroded region are governed by different equations.
in which is the mass density of quartz; 2ℎ′ and ′ are the thickness and the density of the electrodes, respectively; , and are the compact forms of elastic, piezoelectric and dielectric constants; and and are the fundamental TS frequencies for infinite unelectroded and electroded plates respectively, as seen in [18].
To deal with the governing equations we used a method similar to that adopted by He et al. [26]. For both the electroded and unelectroded regions, we introduced the following coordinate transformation:  The thickness-shear vibration consists of different modes. The displacement u n 1 (x 1 , x 3 , t) for the nth-order thickness-shear vibration of AT-cut quartz plate is defined by [5]: where n = 1 represents the fundamental thickness-shear mode, and larger values of n represent higher-order overtone modes. According to [15], the items u n 1 representing vibrations of electroded and unelectroded regions are governed by different two-dimensional scalar differential equations.
The scalar equation for the electroded region is shown as: While the scalar equation for the unelectroded region has the similar form as: where M n = c 11 + (c 12 + c 66 )r + 4(rc 66 − c 66 )(rc 22 + c 12 ) in which ρ is the mass density of quartz; 2h and ρ are the thickness and the density of the electrodes, respectively; c pq , e ip and ε ij are the compact forms of elastic, piezoelectric and dielectric constants; and ω ∞ and ω ∞ are the fundamental TS frequencies for infinite unelectroded and electroded plates respectively, as seen in [18].
To deal with the governing equations we used a method similar to that adopted by He et al. [26]. For both the electroded and unelectroded regions, we introduced the following coordinate transformation: where λ and µ are defined as: and Equations (2) and (3) separately become: The boundary of the electrode can be defined by: After the coordinate transformation the circular boundary becomes an ellipse shown as: and we introduce the elliptical coordinate in order to deal with the boundary conditions: where c is half the focal distance of the ellipse which represents the electrode boundary in elliptical coordinates; and ξ is the radial variable; while η is the angular variable. Then Equations (10) and (11) are transformed into Equations (15) and (16) respectively: in which: We now perform a variable separation to u n 1 : Equations (15) and (16) can be rewritten as: where i = 1, 2; and λ is the separation constant which is related to q. The first and second equations of Equation (19) are known as the Mathieu equation and modified Mathieu equation, and the subscript i = 1, 2 represents the equation used for the electroded and unelectroded parts, respectively.

Mathieu and Modified Mathieu Equations
The exact solutions of the Mathieu and modified Mathieu equations can be obtained from [27]. For the Mathieu equation, there are four kinds of periodic solutions with the period 2π, representing angular distributions, as shown by: Here Equation (20) is called the Mathieu function, in which ce 2m and ce 2m+1 represent symmetric solutions, while se 2m+2 and se 2m+1 represent anti-symmetric solutions. The parameters m = 0, 1, 2 · · · represent solution orders. A m k and B m k are coefficients determined by q. For the modified Mathieu equation there are also four kinds of solutions, representing radial distributions, as shown by: where m = 0, 1, 2 · · · also represents the order of the solutions. Equation (21) is called the modified Mathieu function in which Ce 2m and Ce 2m+1 represent symmetric solutions while Se 2m+2 and Se 2m+1 represent anti-symmetric solutions. A m k and B m k are coefficients, the same as those in Equation (20). All the modified Mathieu functions can be expanded with Bessel functions of the first and second kinds shown as: where v 1 = √ q exp(−ξ); v 2 = √ q exp(ξ); and J = 1, 2. A k m and B k m are mentioned in Equation (21).
, which is a Bessel function of the first kind, and Equation (22) becomes the first kind of modified Mathieu function.
, which is a Bessel function of the second kind, and Equation (22) becomes the second kind of modified Mathieu function. Similar to the definition of the Hankel function, there are also the third and fourth kinds of modified Mathieu functions defined as: Mc Different modified Mathieu functions are used in different kinds of solutions representing different vibration modes [27], which we will discuss in the following sections.

Vibration Modes and Displacement Solutions
Because of the features of the Mathieu and modified Mathieu functions, the vibration modes can be divided into symmetric and anti-symmetric modes.
For symmetric modes, the displacement u n 1 in the electroded region can be written as In the unelectroded region, the displacement u n 1 can be written as For the anti-symmetric modes, the displacement for the electroded region can be written as: In the unelectroded region the displacement is: where A m , B m , A m and B m are undetermined constants. In Equations (24) and (26) we chose the first kind of modified Mathieu function, which represents a standing wave because the electroded region is bounded and there could be a displacement in the middle of the region. In Equations (25) and (27), the third kind of modified Mathieu function, which represents a divergent wave, was chosen because the unelectroded region extends to infinity and the waves can propagate outwards.

Boundary Conditions
We take the symmetric modes as examples. The far-field radiation condition of the unelectroded region has been expressed in [4] by: On the interface between the electroded and the unelectroded regions in x 1 -x 3 plane, the displacement and stress continuity conditions can be expressed as: Equation (28) is taken into account first. The substitution of Equation (25) into Equation (28) results in the following equation: According to the characteristics of the third kind of modified Mathieu function, Mc m approaches zero when ξ = ξ ∞ and Equation (30) is satisfied automatically.
With respect to the continuity conditions, we can substitute Equations (24) and (25)

The Collocation Method
In Equations (28) and (29), the displacements are written as a sum of infinite series of terms and we took several of them as the approximate displacements, and the result converged to the accurate value as more items are considered.
Here we adopted the collocation method, i.e., the continuous equations were satisfied at a series of collocation points on the interface contour between the electroded and the unelectroded parts of the structure in the x 1 -x 3 plane.
We now considered the symmetric modes of the infinite plate. Due to the features of the Mathieu functions and the symmetry of the model, only half the boundary in the first and fourth quadrants needed to be examined for the displacement and stress continuation. After choosing collocation points we found linear homogeneous equations set from Equation (31) for the undetermined constants A m and B m . To ensure the existence of the nontrivial solutions A m and B m , the determinant of the coefficient matrix should be zero, which can determine the value of ω.
In order to use the collocation method, the convergence of the vibration frequency should be first examined to ensure validity.
For example, an infinite plate with an electrode attached on each of the two major surfaces (in the x 1 -x 3 plane) has a thickness of 2h = 0.1 mm, the electrode has a radius of R 0 = 0.8 mm, and the mass ratio between the electrode and the plate is taken to be R = 2ρ h /ρh = 0.05. The material constants of AT-cut quartz can be obtained from [28]. When solving the Mathieu equation we took the first thirty terms of each series to conduct the calculation. Two to nine points on half the boundary were chosen to calculate the frequency ω of the fundamental mode when n = 1, and the result is shown in Figure 2: With respect to the continuity conditions, we can substitute Equations (24) and (25) into Equation (29) to get:

The Collocation Method
In Equations (28) and (29), the displacements are written as a sum of infinite series of terms and we took several of them as the approximate displacements, and the result converged to the accurate value as more items are considered.
Here we adopted the collocation method, i.e., the continuous equations were satisfied at a series of collocation points on the interface contour between the electroded and the unelectroded parts of the structure in the 1 x -3 x plane.
We now considered the symmetric modes of the infinite plate. Due to the features of the Mathieu functions and the symmetry of the model, only half the boundary in the first and fourth quadrants needed to be examined for the displacement and stress continuation. After choosing collocation points we found linear homogeneous equations set from Equation (31)  In order to use the collocation method, the convergence of the vibration frequency should be first examined to ensure validity.
For example, an infinite plate with an electrode attached on each of the two major surfaces (in the 1 x -3 x plane) has a thickness of 2h = 0.1 mm, the electrode has a radius of R0 = 0.8 mm, and the mass ratio between the electrode and the plate is taken to be The material constants of AT-cut quartz can be obtained from [28]. When solving the Mathieu equation we took the first thirty terms of each series to conduct the calculation. Two to nine points on half the boundary were chosen to calculate the frequency  of the fundamental mode when 1 n  , and the result is shown in Figure 2: It can be seen from Figure 2 that when the number of the chosen points increased to five, the frequency calculated became steady and convergence was ensured.
We then checked the displacement and stress continuity of the interface between the electroded and the unelectroded parts on the elliptical coordinate . We arbitrarily chose boundary points between the collocation points for examination. For Example 1, with eight collocation points chosen on half of the interface, and the first eight terms of the series for the calculation of the displacement It can be seen from Figure 2 that when the number of the chosen points increased to five, the frequency calculated became steady and convergence was ensured.
We then checked the displacement and stress continuity of the interface between the electroded and the unelectroded parts on the elliptical coordinate ξ = ξ 0 . We arbitrarily chose boundary points between the collocation points for examination. For Example 1, with eight collocation points chosen on half of the interface, and the first eight terms of the series for the calculation of the displacement and stress remaining, we obtained the frequency of the fundamental TS mode and the expression of the displacement and the stress. The displacement and stress of two arbitrarily-chosen points on the interface between the collocation points (e.g., (ξ 0 ,π/4) and (ξ 0 ,5π/6)) are shown in Table 1. Comparing the data in Table 1, we found that the displacement was continuous at the arbitrarily-chosen points besides the collocated points on the boundary, and the corresponding stress difference across the boundary was also acceptable. Thus, we concluded that the collocation method can be used to deal with the boundary conditions under question.

Numerical Results and Discussion
In this section, we will discuss the displacement distribution of the TS vibrations for Example 1. The displacement field for the fundamental TS mode, u 1 1 (n = 1), is shown in Figure 3a, while the displacement field for the 3rd-order overtone TS mode, u 3 1 (n = 3), is shown in Figure 3b. Comparing Figure 3a,b, we found that the shapes of the distributions were almost the same, but that the displacement u n 1 . was more concentrated in the electroded region when n = 3 than when n = 1. Due to the material anisotropy of the quartz crystal plate, the concentration area of the vibration in the middle of the plate was an ellipse instead of a circle, where the long axis lies in the x 1 direction and the short axis in the x 3 direction. Meanwhile, the frequency for the third thickness-shear overtone was about triple that of the fundamental mode. and stress remaining, we obtained the frequency of the fundamental TS mode and the expression of the displacement and the stress. The displacement and stress of two arbitrarily-chosen points on the interface between the collocation points (e.g., and ) are shown in Table 1. Comparing the data in Table 1, we found that the displacement was continuous at the arbitrarily-chosen points besides the collocated points on the boundary, and the corresponding stress difference across the boundary was also acceptable. Thus, we concluded that the collocation method can be used to deal with the boundary conditions under question.

Numerical Results and Discussion
In this section, we will discuss the displacement distribution of the TS vibrations for Example 1. The displacement field for the fundamental TS mode, 1 1 u (n = 1), is shown in Figure 3a, while the displacement field for the 3rd-order overtone TS mode, 3 1 u ( 3 n  ), is shown in Figure 3b.
Comparing Figure 3a,b, we found that the shapes of the distributions were almost the same, but that the displacement 1 n u was more concentrated in the electroded region when  Four other thickness-shear modes with n = 1 are shown in Figure 4. Among these, Figure 4a,b are modes symmetric to axis- 3 x , while Figure 4c,d are anti-symmetric modes. Figure 4a-d are not the fundamental modes for n = 1 and they have more nodal lines under the electrodes, which will cause charge cancellation and are unwanted spurious modes. In fact, the mode in Figure 3a is the most useful one. Furthermore, by observing the vibration in the electroded and unelectroded areas, it was obvious that the vibration amplitude was low in the unelectroded area and most of the vibration Four other thickness-shear modes with n = 1 are shown in Figure 4. Among these, Figure 4a,b are modes symmetric to axis-x 3 , while Figure 4c,d are anti-symmetric modes. Figure 4a-d are not the fundamental modes for n = 1 and they have more nodal lines under the electrodes, which will cause charge cancellation and are unwanted spurious modes. In fact, the mode in Figure 3a is the most useful one.
Furthermore, by observing the vibration in the electroded and unelectroded areas, it was obvious that the vibration amplitude was low in the unelectroded area and most of the vibration energy was confined to the electroded area. This phenomenon is called energy trapping and it is used for device packaging.

Model, Equation and Boundary Conditions
We then used the COMSOL Multiphysics software to solve the problem in order to verify the numerical results solved by the Mathieu functions in Section 3.
Since the plate is unbounded, we needed to use the infinite element domain whose geometry is set cylindrical. We created an infinite two-dimensional model for the crystal plate in COMSOL. The electroded zone was within the circle at radius 0 R = 0.8 mm, and the outer zone was the unelectroded zone. The outer side of the unelectroded zone tended to infinity. In addition, the two circles outside the electrode were introduced for mesh, the radius of which needed to be big enough to make the model accurate, yet small enough to be economic in computation. After examination, we determined the appropriate size for the outer radiuses of the two rings to be 1.8 mm and 2.5 mm. We used extra fine free triangles in the model to ensure that the results were accurate enough; the meshed model in the 1 x − 3 x plane is shown in Figure 5.
The problem in this article was mode analysis by two-dimensional equations, which is an eigenvalue problem that can be solved by a user-defined Partial Differential Equation (PDE) interface. Coefficient form PDE, which provides a general interface for specifying and solving many well-known PDEs in the coefficient form, was chosen from the PDE interfaces in module mathematics shown as:

Model, Equation and Boundary Conditions
We then used the COMSOL Multiphysics software to solve the problem in order to verify the numerical results solved by the Mathieu functions in Section 3.
Since the plate is unbounded, we needed to use the infinite element domain whose geometry is set cylindrical. We created an infinite two-dimensional model for the crystal plate in COMSOL. The electroded zone was within the circle at radius R 0 = 0.8 mm, and the outer zone was the unelectroded zone. The outer side of the unelectroded zone tended to infinity.
In addition, the two circles outside the electrode were introduced for mesh, the radius of which needed to be big enough to make the model accurate, yet small enough to be economic in computation. After examination, we determined the appropriate size for the outer radiuses of the two rings to be 1.8 mm and 2.5 mm. We used extra fine free triangles in the model to ensure that the results were accurate enough; the meshed model in the x 1 -x 3 plane is shown in Figure 5.
The problem in this article was mode analysis by two-dimensional equations, which is an eigenvalue problem that can be solved by a user-defined Partial Differential Equation (PDE) interface. Coefficient form PDE, which provides a general interface for specifying and solving many well-known PDEs in the coefficient form, was chosen from the PDE interfaces in module mathematics shown as: where u is the unknown scalar to be solved; t and f represent the time and source term, respectively; λ is the eigenvalue; and e a , d a , c, α, γ, β and a are coefficients defined by users in which c is a matrix because of the material anisotropy. This is for a 2D domain, ∇ = [∂/∂x, ∂/∂y].  The results of the fundamental and third overtone TS modes obtained through FEM simulation are shown in Figure 6. Meanwhile, two spurious modes are shown in Figure 7.  The mode frequencies obtained from the theoretical analysis and FEM simulation are also listed in Table 2, where the first and last rows show the fundamental and third order overtone TS modes, respectively, and the four rows in the middle illustrate four typical spurious modes whose frequencies are close to that of fundamental mode ( 1 n  ). We can see that for each mode the error of frequency is less than 100ppm, which is acceptable in the analysis and design of resonators. That validates the correctness of our theoretical analysis based on the coordinate transform and the Mathieu function series expansion.  The results of the fundamental and third overtone TS modes obtained through FEM simulation are shown in Figure 6. Meanwhile, two spurious modes are shown in Figure 7.  The mode frequencies obtained from the theoretical analysis and FEM simulation are also listed in Table 2, where the first and last rows show the fundamental and third order overtone TS modes, respectively, and the four rows in the middle illustrate four typical spurious modes whose frequencies are close to that of fundamental mode ( 1 n  ). We can see that for each mode the error of frequency is less than 100ppm, which is acceptable in the analysis and design of resonators. That validates the correctness of our theoretical analysis based on the coordinate transform and the Mathieu function series expansion.  The mode frequencies obtained from the theoretical analysis and FEM simulation are also listed in Table 2, where the first and last rows show the fundamental and third order overtone TS modes, respectively, and the four rows in the middle illustrate four typical spurious modes whose frequencies are close to that of fundamental mode (n = 1). We can see that for each mode the error of frequency is less than 100 ppm, which is acceptable in the analysis and design of resonators. That validates the correctness of our theoretical analysis based on the coordinate transform and the Mathieu function series expansion.

Conclusions
The exact frequencies of thickness-shear and thickness-twist modes along with mode shapes for infinite, centrally partially-electroded circular AT-cut quartz plates were obtained through theoretical analysis, by the application of coordinate transform and using Mathieu and modified Mathieu functions. The results were compared with 2D FEM simulations by COMSOL software and good agreement was achieved, illustrating the correctness of our theoretical method. From the mode shapes obtained it was easy to find that all the vibrations were concentrated in the electroded region and there was little vibration in the outer, unelectroded area. This phenomenon is called energy trapping and is used for resonator packaging, which does not affect TS vibrations. Our results are of theoretical importance for the design of circular quartz resonators.