Abstract
A general and powerful finite element code has been developed for the analysis of quartz crystal strip resonators utilizing a three dimensional hexahedral element with cubic interpolations of the geometry and field quantities. This higher order interpolation has been found to be superior for computing high frequency vibration modes than the twenty node quadratic elements normally used in three dimensional analysis. The difficulty in designing strip resonators comes mainly from the existence of coupling of unwanted modes such as flexure, face shear, and extension, with the thickness shear mode, and moreover, this coupling is very strongly influenced by crystal blank dimensions. Therefore, in order to properly design these devices, an analysis tool is necessary which can predict these modes with adequate precision. The general development of the theory and finite element equations are presented along with results. An approximate method is employed to remove the electric variables from the formulation, greatly reducing the storage and CPU requirements of the program. A convergence analysis is also presented which quantifies the precision differences between the quadratic and cubic elements.

Introduction
The problem of true three dimensional modeling of high frequency vibrations in crystal resonators has been, and still is, a very challenging problem in computational mechanics. The advent of very powerful, low priced desktop computers in recent years has made three dimensional analysis accessible to design engineers. For this reason, a finite element analysis system for three dimensional (as well as two dimensional) analysis has been developed for use in designing quartz crystal strip resonators for low cost VCXO applications. A general and powerful finite element code has been developed for the analysis of quartz crystal strip resonators utilizing a three dimensional hexahedral element with cubic interpolations of the geometry and field quantities. This higher order interpolation has been found to be superior for computing high frequency vibration modes than the twenty node quadratic elements normally used in three dimensional analysis. The difficulty in designing strip resonators comes mainly from the existence of coupling of unwanted modes such as flexure, face shear, and extension, with the thickness shear mode, and moreover, this coupling is very strongly influenced by the crystal blank dimensions. Therefore, in order to properly design these devices, an analysis tool is necessary which can predict these modes with adequate precision. It has been shown by Yong, et al [3] that the two dimensional first order plate theories do not accurately predict the unwanted modes, and this same author has proposed using higher order plate theory to improve accuracy. The present development takes a different approach, in that a general three dimensional analysis is used to improve the results. This three dimensional analysis represents an alternative methodology which can be compared with the plate theory methods presently under development elsewhere. The general development of the theory and finite element equations are presented along with results from both the two and three dimensional analyses. Convergence analyses are carried out for quadratic and cubic interpolation schemes. The frequency spectrum for a Z' elongated strip resonator is presented and is compared with the equivalent two dimensional results. A powerful method of three dimensional volumetric visualization for mode identification has also been developed. This technique, known as iso-surface extraction, has been investigated previously by the first author [4], and here has been improved for general use in the analysis of vibration mode shapes. The power of this technique lies in its ability to extract a very sparse amount of information from the three dimensional solution model and display a single image for each displacement component from which all necessary information for mode identification can be deduced. Plots of coupled thickness shear and flexural vibration modes are presented illustrating the use of this technique.

Development of The Approximate Three Dimensional Finite Element Equations for Thickness Field Excitation
The development of the three dimensional finite element equations for a piezoelectric crystal plate with parallel conducting electrodes on a portion of its major surfaces begins by considering the general equations of linear piezoelectricity







and upon combining (1) - (3),



and



In equations (1) - (5), T_ij are the components of mechanical stress, D_i are the components of electric displacement, S_ij are the components of strain, E_i are the components of electric field, c_ijkl are the elastic constants, e_kij are the piezoelectric stress constants, _ij are the dielectric constants, and is the mass density of the material. In these same equations, u_k denotes the components of mechanical displacement and represents the electric potential. The variational, or weak form of equations (1) - (5) for a body occupying a volume V bounded by a surface S is given as



where



is the kinetic energy density and



represents the electric enthalpy density. In equation (6) t_i is the applied traction on the bounding surface and D_in_i represents the charge density on the bounding surface with normal n_i.

The approximate equations for thickness field excitation are obtained from the general equations by neglecting the planar variations of the electric variables which leads to the following form of the electric potential for a plate of thickness 2h in the electroded region



and



in the unelectroded region. In equations (9) and (10) the quantities and represent the values of the mechanical displacements on the top (x₂ = +h) and bottom (x₂ = -h) of the plate, respectively. Using (9) and (10), the variation of the electric enthalpy becomes



in the electroded region, and



in the unelectroded region, with





and



representing effective material constants. The effects of thin electrode plating with thickness 2h' and mass density 'on portions of the two major surfaces of the plate can be incorporated into the formulation by imposing the following traction boundary conditions on the portion of the bounding surface containing the electrodes



In equation (16), is the angular frequency and represents the portion of the outer surface, S, containing the electrode plating with the "+" denoting the top surface (x₂ = +h) and the "-" denoting the lower surface (x₂ = -h). Using (11), (12), and (16) in equation (6) gives the final form of the variational equation



The finite element discretization process is applied by interpolating the displacements with a set of shape functions, N^q. as follows



where are the nodal displacements. In equation (18) the superscripts are intended to imply a sum over nodes within a single element or an entire mesh, depending upon context. This notation will be employed throughout to save space. The shape functions N^q. may take on several forms and will not be explicitly defined here. The reader is referred to [1] and [2] for these and other omitted finite element definitions. Using (18) in the variational equation (17) leads to the equation



where



is the usual stiffness matrix found in linear elastic analysis with the new effective elastic constants, and



represents a correction to the stiffness matrix to account for the shorting of the electrodes on the major surfaces of the plate. Likewise,



represents the usual mass matrix, and



is the correction to the mass matrix caused by the addition of thin electrode plating on the major surfaces of the plate. It is worth noting that these two correction matrices are represented by surface integrals rather than the volume integrals used to integrate the normal stiffness and mass matrices. It is also instructive to note that the correction matrices vanish in the unelectroded regions of the plate. Equation (19) can be written in matrix form as



which represents a symmetric, positive definite eigen value problem. The eigen values which solve this system yield the frequencies, and the corresponding eigen vectors determine the mode shapes in the plate associated with these frequencies.

At this point, it is important to appreciate that in the resulting equations (19), the electric variables do not appear. This fact allows for a substantial savings in computational effort at the expense of accuracy. It is not claimed that the approximate method applied here is a substitute for full piezoelectric analysis, but rather a reasonable approximation which greatly simplifies the problem, while maintaining reasonable accuracy for materials with weak electro-mechanical coupling, such as quartz. Furthermore, it must be noted that the present formulation is very specific to the case of a single pair of perfectly parallel and aligned electrodes with equal areas on top and bottom deposited on a flat plate.

Visualization Methods
An integral part of any finite element package is the post processor. Once a solution is obtained, an efficient method is needed to analyze the results. The problems considered here are of such a high frequency that a huge number of eigen pairs are generated in a small frequency interval. This is so because every harmonic of every lower frequency mode is present in the spectrum. In a real crystal plate, the higher harmonics of the lower frequency modes tend to decrease in activity as the harmonic number increases. In the eigen analysis applied here, they do not. What is therefore needed is a method to inspect the vibration modes and determine which ones are important for practical analysis of the resonator. To this end, a three dimensional post processor has been developed which is based on a method of isosurface extraction. Such methods are widely used in medical imaging as well as other areas of scientific visualization. The software that has been developed takes as its input the finite element mesh and solution vector, from which a collection of triangles are created representing the intersection of a surface of constant displacement with the model volume. The triangles can then be quickly rendered by the graphics hardware. The surfaces are extracted from a normalized set of displacements at a threshold between 0.0 and 1.0. For proper visualization of standing wave phenomena, two surfaces are actually extracted, one for the positive value and one for the negative value of the threshold. These surfaces are usually color coded to aid in identification. The advantage of this method over others is that one single image can be used to simultaneously discern vibration characteristics in all three directions.

Analysis of Results
Figure (1) shows a basic schematic layout of a strip resonator. A rectangular crystal plate is mounted on an effectively rigid shelf cantilever style with a single pair of parallel electrodes deposited on the major surfaces. All problems considered in this study correspond to plates of this type.

Shown in figure (2) are the results of a convergence study performed on a 13 Mhz AT Cut strip using two dimensional analysis. In the two dimensional analysis, the cross section of the plate in the x₁ - x₂ plane cut through the center of the electrode is modeled with two dimensional elements. This analysis solves the special case of straight crested waves and, while approximate, is sufficiently accurate to asses the convergence of the calculated frequencies with respect to mesh fineness in the thickness and lateral directions. In this figure the frequencies of a particular strip resonator are plotted for different mesh densities. The x-axis represents the number of elements along the crystal width while the y-axis represents the normalized frequency. Each separate curve represents a different number of elements in the thickness direction. The dashed lines represent quadratic elements using from 2 through 6 elements through the thickness. The solid line represents a mesh of cubic elements using 2 elements through the thickness. It is observed that at approximately 70 elements along x₁ the curve for the cubic elements intersects the curve for 5 quadratic elements. This fineness represents a convergence of better than 1 Khz for this particular plate. Based on these results, it can be concluded that the use of cubic elements is roughly 3 times more accurate than quadratic elements for the convergence in the thickness direction of the plate. This accuracy increase translates into substantial memory savings. For example, at the intersection point noted here, the quadratic mesh contains 917 nodes and 350 elements which results in a storage requirement for the finite element program of 426,749 words. At the same point, the cubic element mesh results in 1201 nodes and 140 elements, which requires 278,365 words of storage for the finite element program, about 65% of the storage required by the quadratic mesh at the same level of accuracy in frequency convergence. In general it has been found that for kilohertz accuracy, approximately two cubic elements per half wave for any displacement component are needed. Generally speaking, the flexural harmonic will govern this convergence in the lateral directions, and the thickness harmonic will govern the convergence in the thickness direction.

Shown in figure (3) is a frequency spectrum for a 10 Mhz AT Cut ( = 35°10') strip resonator. The plate under consideration for this study has the following dimensions (referring to figure (?)): 2l = 5.588 mm, 2h = 0.16346 mm, 2a = 1.7272 mm, 2b = 2.032 mm, d₁ = 1.524 mm, d₂ = 0.762 mm, 2h' = 2804 angstroms gold, and 2w = 2.12498 mm - 2.45190 mm. Plotted on the x-axis of this graph is the aspect ratio of the plate and on the y-axis is plotted the absolute frequency. Shown by the dotted curves are the results of a two dimensional analysis performed on the cross section of the plate, and shown by the solid curves are the results of a three dimensional study done on the entire plate. In the three dimensional analysis, only the family of modes representing transverse anharmonics of even lateral anharmonics of thickness shear coupled with odd lateral anharmonics of flexure are shown. The modes are identified as TS for thickness shear and F for flexure, with the first number denoting the lateral anharmonic and the second number denoting the transverse anharmonic. The first four transverse anharmonics are shown.

In the process of solving for the modes depicted in figure (3), a huge number of eigen pairs were obtained. In order to identify the modes of interest, three dimensional visualization of the mode shapes was necessary. Examples of these plots are shown in figures (4) - (11). Shown in figures (4) - (7) are the u₁ and u₂ displacement components for the first four transverse anharmonics of thickness shear taken at an aspect ratio of l/h = 13.60. In all cases the isosurfaces were thresholded at values of +/-0.4 for u₁ and +/-0.1 for u₂. The "+" and "-" thresholds are always plotted together at different colors to facilitate identification of phase reversals in the mode shapes. Figures (8) - (11) show plots of the same modes at an aspect ratio of l/h = 14.40. Comparing these two sets of modes it is observed that the flexural components in the second group are substantially larger than those in the first. This analysis illustrates the strong dependence of the mode coupling on the width of these plates. The visualization method applied here also allows for the analysis of energy trapping. It can be seen from the plots that the TS1,1 and TS1,2 modes are well trapped along x₃, while the higher modes, TS1,3 and TS1,4 are not, as evidenced by the high displacement activity at the free end of the crystal. The mode shape distribution of the fundamental thickness shear mode along the axis is also sensitive to the crystal width, as shown in figure (12) which plots the mode distributions on the top and bottom surfaces of the plate for different aspect ratios.

Conclusion
The three dimensional finite element code developed here can successfully model the complex behavior of strip resonators. For proper accuracy of the solutions, it is absolutely necessary to use higher order interpolation functions. The cubic elements used here greatly improve the accuracy of the solutions while reducing the memory requirements of the program. It may be so that even higher order elements, such as quintic, may provide even better accuracy. The approximations made to eliminate the electric variables from the formulation reduce the storage requirements of the program by at least 25%. This approximation represents a reasonable trade off with respect to accuracy in the solution.

References
[1] T.J.R. Hughes, The Finite Element Method, Linear Static and Dynamic Analysis, Prentice Hall, 1987.

[2] O.C. Zienkiewicz and R.L. Taylor, The Finite Element Method, Vols. 1&2, McGraw-Hill, 1989.

[3] Y.K. Yong, Z. Zhang, and J. Hou, "On The Accuracy of Plate Theories for the Prediction of Unwanted Modes Near the Fundamental Thickness Shear Mode, Proceedings of the 1995 IEEE Frequency Control Symposium, pp. 755-760.

[4] "Visualization Techniques for Bulk Wave Resonators", D. Silver, J.T. Stewart, Y.K. Yong, Proc. 1991 IEEE Ultrasonics Symposium, pp 499-504.

Figures