Resources

Three Dimensional Finite Element Modeling of Quartz Crystal Strip Resonators

J.T. Stewart and D.S. Stevens

Vectron International, 267 Lowell Road, Hudson, NH 03051

**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_{i}n_{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_{2} = +h) and bottom (x_{2} = -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 _{2} = +h*)
and the "-" denoting the lower surface
(

The finite element discretization process is applied by interpolating the displacements with a set of shape functions, N

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

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.

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.

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

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,

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

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.

[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.

Figure 1. Basic schematic layout of a crystal strip resonator.

Figure 2. Convergence analysis of quadratic and cubic elements.

Figure 3. Frequency spectrum for a 10 Mhz AT-Cut quartz strip resonator.

Figure 4. TS1,1 mode l/h = 13.60.

Figure 5. TS1,2 mode l/h = 13.60.

Figure 6. TS1,3 mode l/h = 13.60.

Figure 7. TS1,4 mode l/h = 13.60.

Figure 8. TS1,1 mode l/h = 14.40.

Figure 9. TS1,2 mode l/h = 14.40.

Figure 10. TS1,3 mode l/h = 14.40.

Figure 11. TS1,4 mode l/h = 14.40.

Figure 12.