Helping Customers Innovate, Improve & Grow

J.T. Stewart and D.S. Stevens

Vectron International, 267 Lowell Road, Hudson, NH 03051

**Abstract**

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.

**Introduction**

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 [4] 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

where

and

.

In equations (1)-(4)
T_{ij}
are the components of the stress tensor,
s_{ij}
are the components of the infinitesimal strain tensor,
u_{j}
are the components of displacement,
b_{j}
are the components of the body force per unit volume,
c_{ijkl}
are the components of the elastic stiffness tensor, and
v_{ij}
are the thermoelastic constants with
_{ij}
representing the thermal expansion constants for the
material. In equation (3) the quantity
(T - T_{0})
represents the small variation of temperature from the ambient value,
T_{0}.

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,
u_{i}
are defined in the usual way and
t_{i}
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,
N^{q}
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
N^{q}
may take on several forms
and will not be explicitly defined here. The reader is
referred to [10] and [11] 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

where

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

with solution

**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
[3]. 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
x_{1}
length 2a,
x_{3}
width 2b, and thickness 2h as

and

The various constants appearing in equations (13) - (24)
are given as

In the above, the quantities
^{(1)},
^{(2)},
^{(3)},
r_{2},
r_{3},
_{2},
and

_{3},
are defined in [3]. In addition, the constants
c_{16},
c_{12},
c_{57},
appearing in equations (25) - (31) are not
the ordinary elastic constants but rather special derived
constants which are also defined in [3]. The reader is
referred to [3] for a detailed discussion on the derivation of
all of these quantities. The amplitude constants appearing
in equations (13)-(24) are given as

and

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 [3]. As discussed in [3], 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
x_{1} - x_{3}
plane, defined by the components
R_{ij},
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,
(x_{1},x_{2},x_{3}),
to a new system
(x'_{1},x_{2},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 [2] for small fields superposed on a bias [1].
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

where

and

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 - T_{0}),
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 [7]. 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
n^{th}
element. A single
element integral as defined by equation (46) is evaluated by
first subdividing the element domain,
_{n},
into an
n_{1} x n_{2} x n_{3}
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
(i^{th}, j^{th}, k^{th})
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
M_{1} x M_{2} x M_{3}
Gaussian Quadrature rule

W^{m}
1
W^{m}
2
W^{m}
3

wherein the mode shape derivatives,
,
are evaluated at each
Gauss point

(),
and summed along with the proper weights,
W^{m}
,
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,
*t _{f}* = 3.0 mm,

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.

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

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

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

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

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

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

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

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

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

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

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

Figure 1. Four regions for the definition of the trapped energy mode shape.

Figure 2. Angle definition for support rotation studies.

Figure 3. Schematic diagram for shelf type mount.

Figure 4. Finite element model for shelf type mount.

Figure 5. Sensitivity vs. support orientation for shelf type mount.

Figure 6. Temperature compensated angle vs. orientation for shelf type mount.

Figure 7. Schematic diagram for edge type mount.

Figure 8. Finite element model for edge type mount.

Figure 9. Sensitivity vs. support orientation for edge type mount.

Figure 10. Sensitivity vs. support orientation for edge type mount for higher overtones.

Figure 11. Schematic diagram for clip mount.

Figure 12. Finite element model for clip mount.

Figure 13. Sensitivity vs. support orientation for clip mount for various clip lengths.

Figure 14. Schematic diagram for shelf type mount with inverted mesa crystal.

Figure 15. Finite element model for shelf type mount with inverted mesa crystal.

Figure 16. Comparison of inverted mesa and flat crystals with shelf type mount.