A general finite element method has been developed for the analysis of the effects of thermally induced mounting stresses on the resonant frequency of crystal resonators. The present study analyses this behavior using a finite element model to obtain an accurate description of the thermal stresses and deformations in a crystal plate that has been mounted in various configurations using an adhesive to bond the quartz to a ceramic base. This solution is then combined with an analytical solution of the vibration mode shape for a harmonic of thickness shear in a general perturbation procedure to obtain the frequency shift. The methods developed herein are useful for studying the sensitivity of a particular mounting scheme to thermal stress for small temperature excursions. Examples of circular AT Cut crystal plates mounted in various configurations are studied.
When a crystal resonator is mounted in a package, stresses are imparted to it from the fixity method which has been employed, which cause the frequency to shift by a small amount. These stresses can be immediate, or can evolve slowly over time, giving rise to a portion of the aging characteristics observed in the device. Mounting stresses affect the frequency of a crystal resonator in two basic ways. The first effect is manifested when the crystal is mounted to another structure with differing thermal expansion properties, causing thermal stresses to develop when the temperature is changed. The second effect comes about when a crystal blank is bonded to a structure with an adhesive at a temperature different than the normal operating temperature of the device, causing a residual stress to develop at the ambient temperature. Depending on the bonding agents used, this residual stress will then "relax" over time, giving rise to an aging effect in the device. The present study seeks to analyze the first of these effects using a finite element model to obtain an accurate description of the thermally induced stresses and deformations in a crystal plate that has been mounted in various configurations using an adhesive to bond the quartz to a ceramic base. This solution is then combined with an analytical solution of the vibration mode shape for a harmonic of thickness shear in a general perturbation procedure to obtain the frequency shift. The perturbation procedure used here employs the first order temperature coefficients of the crystal in addition to the effective constants resulting from the stress bias. This method is an approximation to the more general higher order case which is accurate for small temperature excursions.
A finite element model has been created which allows the properties of the epoxy bond to be included, and a general analytical solution of a trapped energy thickness shear mode for a doubly rotated quartz crystal plate has been developed. A general numerical integration procedure is used to combine an arbitrary mode shape with any static finite element solution to obtain the first perturbation of the eigen-value as outlined by Tiersten's theory [1,2]. This program is used to analyze the mounting stress effects discussed above as a function of crystal cut and mounting configuration for circular AT cut crystals. To compute the frequency shift, a general purpose integration method has been developed which calculates the first perturbation of the eigen-value from a finite element solution of the static biasing state and an analytical solution of the mode shape for a trapped energy resonator operating in an arbitrary overtone of thickness shear. The numerical technique that has been developed generalizes the methods employed by L.D.Clayton and E.P. Eernisse  with many enhancements to improve accuracy. The mode shape that has been developed is completely general accommodating either singly or doubly rotated cuts of quartz through the well know planar transformation of the mode shape axes. The integration procedure readily accommodates this planar transformation by using a very general quadrature procedure combined with a sub-meshing algorithm which decomposes an arbitrarily defined finite element (linear, quadratic, cubic,...) into a set of linear hexahedral integration cells in which the integral is evaluated. Previous work by the first author in SAW acceleration sensitivity [5,6] used a simpler scheme which divides the domain into cube-shaped cells and evaluates the desired quantity exactly in each sub-region. This method is not generally applicable to BAW problems because of the possible planar rotation of the mode-shape axes and the existence of non-rectangular geometries caused by circular or elliptical electrodes and/or contouring. The present method overcomes these difficulties by using a sophisticated quadrature method for evaluating the integral of the product of the mode shape derivatives in a general linear hexahedron.
Solution of the Biasing State
The development of the static thermoelastic finite element equations for the solution of the biasing state begins with the general three dimensional equations of thermoelasticity
In equations (1)-(4) Tij are the components of the stress tensor, sij are the components of the infinitesimal strain tensor, uj are the components of displacement, bj are the components of the body force per unit volume, cijkl are the components of the elastic stiffness tensor, and vij are the thermoelastic constants with ij representing the thermal expansion constants for the material. In equation (3) the quantity (T - T0) represents the small variation of temperature from the ambient value, T0.
The variational or weak form of equations (1)-(4) is formulated for a body occupying a volume V bounded by a surface S as
where the variational displacements, ui are defined in the usual way and ti represents prescribed tractions on the surface of the body. The finite element discretization process is applied by interpolating the displacements with a set of shape functions, Nq as follows
where are the nodal displacements. In equation (6) 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 Nq may take on several forms and will not be explicitly defined here. The reader is referred to  and  for these and other omitted finite element definitions. Using equation (6) in the functional (5) gives, in the absence of body forces and applied surface tractions,
Here, represents the discretized domain with bounding surface . For arbitrary variations , equation (7) reduces to
is the usual stiffness matrix, and
represents the effective nodal force vector. This force vector represents a set of loads applied to each node in the mesh to produce the strains which are compatible with the thermal expansion of the material. The global finite element matrix system is assembled in the usual way giving rise to the matrix problem
Solution of the Mode Shape
The analytical solution of the trapped energy mode shape in a flat plate with rectangular electrodes is derived from the basic solution by Stevens and Tiersten . This solution represents a simple approximate analytical expression for the displacement field of an arbitrary harmonic of thickness shear with transverse variation. The solution is defined piecewise in four basic regions, denoted by "-", "s", "T", and "c", as depicted in figure (1). The components of displacement in each region are given in general terms for a rectangular electrode of x1 length 2a, x3 width 2b, and thickness 2h as
The various constants appearing in equations (13) - (24) are given as
In the above, the quantities (1), (2), (3), r2, r3, 2, and
3, are defined in . In addition, the constants c16, c12, c57, appearing in equations (25) - (31) are not the ordinary elastic constants but rather special derived constants which are also defined in . The reader is referred to  for a detailed discussion on the derivation of all of these quantities. The amplitude constants appearing in equations (13)-(24) are given as
with an arbitrary constant, which for most purposes can be taken as unity. The quantities , , , and , are propagation and decay constants for a particular trapped energy anharmonic which are determined numerically by the solution of a series of transcendental equations. The general procedure for obtaining these quantities can be found in . As discussed in , the displacements, as given by equations (13) - (24), are represented in a different coordinate system than the axes of the model. These displacements are in fact referred to the thickness solution eigen vector system via a transformation, ir. The relationship between the original coordinate system and the transformed system is given as
Furthermore, there may be a planar rotation of the mode shape in the x1 - x3 plane, defined by the components Rij, which also has to be taken into account. In general, the solution of a given anharmonic of thickness shear will have a mode distribution which is transformed from the original coordinate system, (x1,x2,x3), to a new system (x'1,x2,x'3), rotated relative to the original coordinate system in which the finite element model is represented. The relationship between the rotated and unrotated systems is given as
In summary, to obtain the original displacements referred to the same coordinate system as the finite element model, the two transformations (37) and (36) must be employed. The final step in the calculation of the mode shape, is the mode normalization required by the perturbation process which is given as
where the normalization constant is given as
with the mass density of the material.
Calculation of The Frequency Shift
The frequency shift under the action of a given static biasing state is computed using Tiersten's perturbation method  for small fields superposed on a bias . The procedure used here follows closely the methods employed in [8,9] which considered problems of thermally induced stresses caused by thick electrodes. In general, the change in resonant frequency, , of the th, eigen mode, at frequency , is given as
In equation (41), are the spatially varying effective elastic constants derived from the biasing state with the biasing stresses, the biasing strains, and the biasing deformation gradients. The components represent the effective elastic constants at the temperature T as defined by equation (45). This definition of is valid for small temperature deviations. In general, these constants are a nonlinear function of (T - T0), and for large temperature excursions, the second and third order temperature derivatives are required, in addition to the first order term, , contained in equation (45). The values for used in the present study were obtained from . At this point it should also be mentioned that equation (41) is valid only for small temperature excursions. For larger temperature excursions in the presence of thermal stresses, higher order elastic constants as well as the temperature derivatives of the third order elastic constants are necessary. In equation (42), are the third order elastic constants for the material and are the thermoelastic constants defined by equation (4). In performing the mode shape derivatives, , appearing in equation (41), care must be taken to properly represent these quantities in light of the mode shape rotation defined by equation (37).
The perturbation integral (41) is evaluated as a sum over N elements in the mesh as
where n denotes the volumetric domain of the nth element. A single element integral as defined by equation (46) is evaluated by first subdividing the element domain, n, into an n1 x n2 x n3 grid of linear hexahedric integration cells and sampling the spatially varying components of at the centroid of each cell. For a sufficiently small cell, the values of can be assumed to be approximately constant and can thus be removed from the subintegral. Using this assumption, each term in the sum of equation (46) can be written as
where ((i), (j), (k)) are the coordinates of the centroid of the (ith, jth, kth) integration cell, n(i,j,k). The integral of the normalized mode shape derivatives appearing on the right hand side of equation (47) is evaluated inside the linear hexahedral integration cell using the M1 x M2 x M3 Gaussian Quadrature rule
Wm 1 Wm 2 Wm 3
wherein the mode shape derivatives, , are evaluated at each Gauss point
(), and summed along with the proper weights, Wm , to provide a very accurate numerical integration of these quantities. The quantity J(), appearing in equation (48) is the Jacobian of the integration cell evaluated at the respective Gauss point.
Analysis of Results
Using the developed program, four example problems are studied. In all cases, a 10 Mhz circular AT Cut crystal is mounted in various configurations, and the basic sensitivity of the mount to thermal stresses is computed for different bonding orientations. In all the studies presented, the crystal blank is rotated on its mount through various angles as depicted in figure (2). Figure (3) shows schematically the configuration of a circular crystal mounted to a rectangular ceramic frame with a shelf type bond. The dimensions of the problem under consideration are as follows: 2l = 8.89 mm, 2a = 4.445 mm, 2h = 0.1644 mm, tf = 3.0 mm, hf = 2.0 mm, hg = 0.25 mm. Shown in figure (4) is a typical finite element model of this structure. Figure (5) shows the relationship between the calculated frequency shift and support orientation angle for this problem using a crystal cut of = 35°15'. The one point mount case shows almost no variation with an almost zero frequency shift, which is the expected result for this case, since it approximates very well the case of free expansion. The two point mount shows a wide variation in sensitivity with an apparent compensated orientation of approximately = 63°. This result is consistent with the well known force-frequency behavior for this particular cut. The three and four point mounts do not exhibit a compensating orientation. Shown in figure (6) are the results of an analysis of the same problem in which at each orientation, the frequency shift is computed at slightly different crystal cuts, from which the angle is determined that gives a zero frequency shift at temperature T0. The results are presented as a difference from = 35°15' in minutes. Figure (7) shows schematically the configuration of a circular crystal mounted to a rectangular ceramic frame along the edge of the blank. The dimensions of the problem under consideration are as follows: 2l = 8.89 mm, 2a = 4.445 mm, 2h = 0.1644 mm, tf = 3.0 mm, hf = 2.0 mm, hg = 0.25 mm, and tg = 0.25 mm. Shown in figure (8) is a typical finite element model of this structure. Figure (9) shows the results of the analysis performed for this mount which follow very closely the results of the shelf type mount. Figure (10) shows the results of this same analysis performed for first, third, and fifth overtone operation of the crystal. It is observed that the basic functional relationship between the support orientation angle and the frequency shift is preserved with only slight differences noticed for the higher harmonics. Figure (11) shows schematically the configuration of a circular crystal mounted with straight clips. The dimensions of the problem under consideration are as follows: 2l = 8.89 mm, 2a = 4.445 mm, 2h = 0.1644 mm, hg = 0.20 mm, tg = 0.20 mm, and lc = 0.5 mm - 1.5 mm. Shown in figure (12) is a typical finite element model of this structure. Figure (13) shows the results of an analysis performed for this configuration using a two point mount for various clip lengths. It is observed that the compensating orientation at = 63°, remains the same for each length, while the overall sensitivity increases with shorter clips which are less compliant. The last example, as shown in figures (14) and (15), considers a circular crystal with a shelf type mount on a ceramic frame similar to the first example, except that an inverted mesa crystal is used in place of the flat blank. All dimensions in this example are the same as those of the first, with the addition of the mesa height of hm = 0.25 mm. Figure (16) shows the relationship between support orientation angle and frequency shift for the inverted mesa blank compared to the flat blank for a two point mount. It is observed that there is a roughly 10 degree shift in the compensating orientation, as well as a difference in the overall sensitivity between these two crystals.
The developed finite element software represents an excellent modeling tool which can aid in the design of crystal resonators. The finite element method is certainly the method of choice for this work since a large number of different mount geometries can be studied with the same program, as evidenced by the examples shown. The use of the sub-meshing algorithm with the quadrature integration method, greatly improve the accuracy of the results with a high degree of computational efficiency.
 J.C. Baumhauer and H.F. Tiersten, "Nonlinear Electroelastic Equations for Small Fields Superposed on a Bias", J. Acoust. Soc. Am., Vol. 54, No. 4, 1973, pp. 1017-1034.
 H.F. Tiersten, "Perturbation Theory For Linear Electroelastic Equations for Small Fields Superposed On a Bias", J. Acoust. Soc. Am., Vol. 64, No. 3, Sept., 1978 pp. 832-837.
 D.S. Stevens and H.F. Tiersten, "An Analysis of Doubly Rotated Quartz Resonators Utilizing essentially Thickness Modes With Transverse Variation", J. Acoust. Soc. Am., Vol. 79, No. 6, June 1986, pp. 1811-1826.
 L.P. Clayton and E.P. Eernisse, "Frequency Shift Calculations for Quartz Resonators",Proceedings of the 1991 IEEE International Frequency Control Symposium, pp. 309-320.
 J.T. Stewart, R.C. McGowan, J.A. Kosinski, and A. Ballato, "Semi-Analytical Finite Element Analysis of Acceleration-Induced Frequency Change In SAW Resonators", Proc. 1995 IEEE International Frequency Control Symposium, pp. 499-506.
 J.T. Stewart, J.A. Kosinski, and A. Ballato, "An Analysis of The Dynamic Behavior and Acceleration Sensitivity of a SAW Resonators Supported By Flexible Beams", Proc. 1995 IEEE International Frequency Control Symposium, pp. 507-513.
 B.K. Sinha and H.F. Tiersten, "First Temperature Derivatives of the Fundamental Elastic Constants of Quartz", J. Appl. Phys., 50(4) April 1979, pp. 2732-2739.
 H.F. Tiersten and B.K. Sinha, "Temperature Dependence of the Resonant Frequency of Electroded DoublyRotated Quartz Thickness Mode Resonators", J. Appl. Phys., 50(12) December,1979, pp. 8038-8051.
 D.S. Stevens and H.F. Tiersten, "Temperature Dependence of the Resonant Frequency of Electroded Contoured AT-cut Quartz Crystal Resonators", J. Appl. Phys., 54(4) April,1983, pp. 1704-1716.
 T.J.R. Hughes, The Finite Element Method, Linear Static and Dynamic Analysis, Prentice Hall, 1987.
 O.C. Zienkiewicz and R.L. Taylor, The Finite Element Method, Vols. 1&2, McGraw-Hill, 1989.