Helping Customers Innovate, Improve & Grow
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)
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 [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
x1
length 2a,
x3
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),
r2,
r3,
2,
and
3,
are defined in [3]. 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 [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
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 [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 - 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 [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
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.
Conclusion
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.
References
[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.
Figures
