INTERNATIONAL JOURNAL FOR NUMERICAL METHODS IN ENGINEERING

Int. J. Numer. Meth. Engng (2008)

Published online in Wiley InterScience (www.interscience.wiley.com). DOI: 10.1002/nme.2502
A nite element formulation for thermoelastic damping analysis
Enrico Serra
1, ,
and Michele Bonaldi
2, 3
1
Fondazione Bruno Kessler, FBK-irst MicroTechnologies Laboratory, via Sommarive 18,
I-38100 Povo (Trento), Italy
2
Istituto di Fotonica e Nanotecnologie CNR-FBK, via alla Cascata, 56/C 38100 Povo (Trento), Italy
3
Istituto Nazionale di Fisica Nucleare, via Sommarive 14, I-38100 Povo (Trento), Italy
SUMMARY
We present a nite element formulation based on a weak form of the boundary value problem for fully
coupled thermoelasticity. The thermoelastic damping is calculated from the irreversible ow of entropy
due to the thermal uxes that have originated from the volumetric strain variations. Within our weak
formulation we dene a dissipation function that can be integrated over an oscillation period to evaluate the
thermoelastic damping. We show the physical meaning of this dissipation function in the framework
of the well-known Biots variational principle of thermoelasticity. The coupled nite element equations
are derived by considering harmonic small variations of displacement and temperature with respect to the
thermodynamic equilibrium state. In the nite element formulation two elements are considered: the rst
is a new 8-node thermoelastic element based on the ReissnerMindlin plate theory, which can be used for
modeling thin or moderately thick structures, while the second is a standard three-dimensional 20-node
iso-parametric thermoelastic element, which is suitable to model massive structures. For the 8-node
element the dissipation along the plate thickness has been taken into account by introducing a through-
the-thickness dependence of the temperature shape function. With this assumption the unknowns and the
computational effort are minimized. Comparisons with analytical results for thin beams are shown to
illustrate the performances of those coupled-eld elements. Copyright q 2008 John Wiley & Sons, Ltd.
Received 9 May 2008; Revised 7 October 2008; Accepted 9 October 2008
KEY WORDS: nite elements; thermoelastic damping; ReissnerMindlin plate theory
1. INTRODUCTION
Thermoelastic damping is an internal friction mechanism that originates from stress inhomo-
geneities which in turn generate heat uxes that increase the entropy of a vibrating solid. In MEMS

Correspondence to: Enrico Serra, Fondazione Bruno Kessler, FBK-irst MicroTechnologies Laboratory, via Sommarive
18, I-38100 Povo (Trento), Italy.

E-mail: [email protected]
Contract/grant sponsor: European Community; contract/grant number: RII3-CT-2004-506222
Contract/grant sponsor: Provincia Autonoma di Trento
Copyright q 2008 John Wiley & Sons, Ltd.
E. SERRA AND M. BONALDI
and NEMS devices, thermoelastic damping establishes an absolute lower bound on internal friction
and, in many cases, is the dominant component of damping at room temperature (300 K). For
this reason, the study of thermoelastic damping in exures is an active area of current theoretical
research [1, 2]. In precision measurements, thermoelastic damping acts as a source of mechanical
thermal noise and contributes to limit the sensitivity of gravitational wave detectors [3, 4].
The thermoelastic damping was studied rst by Zener [57]. In particular, he considered the
thermoelastic damping in beams and he gave an analytical approximation for the loss angle of
those structures. More recently, Lifshitz and Roukes [1] obtained a closed form of the thermoelastic
damping problem in a thin beam of rectangular cross section under exural vibrations. They
calculated the loss angle of microbeams and found results close to those obtained with Zeners
model. Another simple formula based on the extension of the Lifshitz and Roukes approach was
obtained by Wong et al. [8] for the computation of thermoelastic damping of MEMS ring resonators.
Nayfeh and Younis [9] developed an independent approach for the evaluation of the loss angle of
microplates: they considered the biharmonic differential plate equation for small transverse (out-
of-plane) displacement of the thin plate together with the heat equation. The thermoelastic damping
was then computed by a perturbative analysis starting from the normal modes of the plate.
The ability to predict thermoelastic damping of all these analytical approaches fails when more
complex geometries or boundary conditions are considered. In fact, the loss angles of structural
modes, which vary spatially across the plate width, cannot be predicted using a simple beam or
plate model. In such cases, a numerical procedure may overcome this limitations.
Recently some efforts have been done in developing nite elements for thermoelastic analysis.
The nite element equations corresponding to the generalized thermoelasticity with two relaxation
times were derived and solved directly in time-domain [10] to investigate second sound effect of
heat conduction in solids subject to thermal shock loading. The thermoelastic damping problem
involved in beam resonators has been addressed using a combination of nite element method
and eigenvalue formulation [11]. This approach allows the evaluation of the eigenfrequencies and
their associated quality factor of resonance, even in case of complex geometries and boundary
conditions. However, a complete framework for the nite element computation of the thermoelastic
damping involving plates and/or massive bodies is still lacking and the physical context is not fully
described. The purpose of this paper is to present a new nite element formulation for computing
the loss angle of a vibrating structure using either planar or full three-dimensional elements. We
derive a weak form for linearly coupled thermoelasticity boundary value problem (BVP), where
an innite velocity of the thermal eld and small harmonic temperature variations with respect to
the equilibrium temperature are assumed. From the weak form a dissipation function is recognized
and used for evaluating the work lost due to thermoelastic effect. We show that the dissipation
function is substantially the same of the dissipation function introduced by M. A. Biot in its early
variational principle [12, 13]. On the other hand our choice of variables leads to a nite element
scheme which requires less DOFs with respect to a nite element implementation based on the
M. A. Biot variational principle. In the nite element formulation two elements are considered.
The rst is a new 8-node plane ReissnerMindlin element which can be used in modeling thin and
moderately thick structures like coatings or the interface layer between two bonded materials. The
second thermoelastic element is a 20-node element and can be used to model thick structures. For
the 8-node element we use the well-known ReissnerMindlin assumptions and the temperature
eld along the element thickness is taken into account by properly modifying the element shape
functions. To validate the procedure we compute the thermoelastic damping of several vibrational
modes of a thin beam and show the comparison with the analytical values obtained by Lifshitz
Copyright q 2008 John Wiley & Sons, Ltd. Int. J. Numer. Meth. Engng (2008)
DOI: 10.1002/nme
A FINITE ELEMENT FORMULATION FOR THERMOELASTIC DAMPING
and Roukes [1]. The outline of the paper is as follows. Section 2 gives an overview of the basic
equations for thermoelasticity, presenting a brief summary of the balance and of the constitutive
equations used in the nite element formulation. We introduce here the dissipation function used
to evaluate the thermoelastic damping. In Section 3 we present the weak formulation of coupled
thermoelasticity BVP and show the connections with the well-known Biots variational principle.
Section 4 introduces the thermoelastic elements with their constitutive relations. The general Finite
Element framework is detailed in Section 5, where the numerical procedure used to compute the
thermoelastic loss angle is also illustrated. In Section 6 analytical results are compared with the
numerical ones for the dissipation of a thin beam. The paper ends with concluding remarks in
Section 7.
2. MATHEMATICAL FRAMEWORK FOR THE THERMOELASTIC DAMPING
We now review the BVP for linearly coupled thermoelasticity in case of small vibrations and small
temperature changes with respect to the equilibrium temperature T
0
. Straindisplacement relation,
balance equations and constitutive equations are needed to formulate the general thermoelastic
problem. We also recall the denition of thermoelastic damping of a body.
2.1. Basic equations for linear coupled thermoelasticity in solids
2.1.1. Straindisplacement equation. When considering small strains and displacements, the kine-
matics equations consist of the well-known straindisplacement relation:
c
i j
=
u
i, j
+u
j,i
2
(1)
2.1.2. Balance equations. Any state of the body must satisfy the balance equation of entropy and
the equation of motion.
The local balance equation of entropy is
T
0
s =q
i,i
(2)
where s is the entropy rate per unit volume, T
0
the equilibrium temperature, q
i,i
the divergence of
the thermal ux. We assumed small temperature changes close to the equilibrium temperature T
0
to obtain a linear equation. No internal heat sources are considered in Equation (2).
The equation of motion in local form is
o
i j, j
+ f
i
=j u
i
(3)
where o
i j, j
is the derivative of the stress tensor, f
i
the volume forces per unit volume, j the density
and u
i
the displacement vector.
2.1.3. Constitutive equations. In thermodynamics, the state of a solid is entirely determined by
the values of a certain set of independent variables: the kinematic variables (strains c
i j
, . . .) and
the temperature. Starting from the rst and the second law of thermodynamics and, for instance,
following the approach in [14], it is possible to show that the Cauchy stress tensor o
i j
is conjugated
Copyright q 2008 John Wiley & Sons, Ltd. Int. J. Numer. Meth. Engng (2008)
DOI: 10.1002/nme
E. SERRA AND M. BONALDI
to the strain tensor c
i j
, in the same way the entropy per unit volume s is thermodynamically
conjugated to the temperature shift 0=T T
0
:
o
i j
=
*F
*c
i j
, s =
*F
*0
(4)
where F(c
i j
, 0) is Helmholtzs free energy. If we retain only the quadratic terms in the Taylor
series expansion in the vicinity of the equilibrium state (c
i j
=0, 0=0), we obtain
F(c
i j
, 0) =
1
2
C
i j kl
c
i j
c
kl
[
i j
c
i j
0+
c
E
2T
0
0
2
(5)
where
C
i j kl
=
*
2
F(0, 0)
*c
i j
*c
kl
, [
i j
=
*
2
F(0, 0)
*c
i j
*0
, c
E
=T
0
*
2
F(0, 0)
*0
2
(6)
In deriving (5) we assume that the free energy, entropy and stress vanish in the equilibrium
state: F(0, 0) =0, s(0, 0) =0, o
i j
(0, 0) =0. From (4) and (5) the two constitutive equations can be
obtained
o
i j
=
*F
*c
i j
=C
i j kl
c
kl
[
i j
0 (7)
s =
*F
*0
=[
i j
c
i j
+
c
E
T
0
0 (8)
where c
E
is the specic heat for the unit volume in the absence of deformation.
The constitutive relation between the heat ux and the temperature gradient is given by Fouriers
law of heat conduction:
q
i
=k
i j
0
, j
(9)
where k
i j
is the thermal conductivity tensor.
In the following we limit our study to the case of isotropic body, where:
C
i j kl
= zo
i j
o
kl
+j(o
i k
o
j l
+o
il
o
j k
)
[
i j
= :(3z+2j)o
i j
k
i j
= o
i j
(10)
Here, : is the coefcient of thermal expansion, z, j are the Lam e constants and is the thermal
conductivity.
2.2. Thermoelastic damping
According to the approach developed by Liftshitz and Roukes [1], the effect of thermoelastic
coupling on the vibrations of a thin beam can be evaluated by solving the coupled equation of
motion to obtain the normal modes of vibration and their corresponding eigenfrequencies c
n
. In
general, the eigenfrequencies are complex and their real part gives the oscillation frequency while
Copyright q 2008 John Wiley & Sons, Ltd. Int. J. Numer. Meth. Engng (2008)
DOI: 10.1002/nme
A FINITE ELEMENT FORMULATION FOR THERMOELASTIC DAMPING
the imaginary part gives the attenuation of the vibration. The loss angle [ at the frequency c
n
is
then obtained as:
[=2

Im(c
n
)
Re(c
n
)

(11)
The loss angle [ can be also evaluated as the work lost per radiant normalized to the maximum
elastic energy stored in body during a cycle. In this case
[=
W
2W
m
(12)
where W is the total work lost per cycle and W
m
is the maximum strain energy during that cycle.
The work lost per cycle W is related to the irreversible entropy produced per cycle (GouyStodola
theorem):
W =
_
cycle
_
V
0
,i
k
i j
0
, j
T
0
dVdt (13)
as can be easily shown by combining Equations (2), (9) and after integrating over the bodys
volume [12, 15, 16]. If we dene the dissipation function D as
D=
1
2
_
V
0
,i
k
i j
0
, j
dV (14)
the loss angle [ may be written as:
[=
_
Ddt
T
0
W
m
(15)
2.3. Two-way and one-way thermoelastic coupling
From Equation (2), using the entropy constitutive equation (8) and denitions (10) for an isotropic
body, we derive the classical linear heat differential equation that can be found in many textbooks:
0
,i i
=c
E
*0
*t
+
E:T
0
(12v)
*c
kk
*t
(16)
where c
E
, E, v is the volumetric specic heat per unit volume in the absence of deformation,
Youngs modulus of the body and the Poisson ratio. Small temperature variations are assumed and
the term related to thermoelastic damping contains time variation of the strain trace c
kk
.
Equivalently [17] the source term can be written in terms of the stress trace o
kk
and the heat
equation becomes:
0
,i i
=c
E
*0
*t
+:T
0
*o
kk
*t
(17)
These equations take into account that the strain eld or the stress eld affects the temperature
eld and conversely the temperature eld affects the strain eld or the stress eld. This is known
as two-way coupling. In many cases the thermoelastic damping is evaluated by considering the
strain or the stress trace in isothermal condition. This approximation is called one-way coupling
because the heat equation is solved with a source term (strain or stress trace) which is independent
from the temperature eld [17].
Copyright q 2008 John Wiley & Sons, Ltd. Int. J. Numer. Meth. Engng (2008)
DOI: 10.1002/nme
E. SERRA AND M. BONALDI
3. FINITE ELEMENTS FOR THERMOELASTIC DAMPING ANALYSIS
We now overview the weak form for thermoelasticity and the approach of M. A. Biot, showing
its connections with the dissipation function D.
3.1. Derivation of the weak form for linear coupled thermoelasticity
We derive the matrix equations for the nite element implementation by using a weak form of
the thermoelastic BVP. Consider a solid at thermodynamic equilibrium, described by the state
variables (u
i
, 0) and identied with the closed domain V in Figure 1. Consider now a variation
(ou
i
, o0) of the state variables around their equilibrium value. The local balance equation (2) can
be multiplied by o0 and integrated over the region V of the body, to obtain:
_
V
(T
0
s +q
i,i
)o0dV =0 (18)
On the second term of this equation we apply the divergence theorem:
_
V
q
i,i
o0dV =
_
V
q
i
o0
,i
dV +
_
A
q
i
n
i
o0dA (19)
then the integral equation (18) may be rewritten as:
_
V
T
0
so0dV
_
V
q
i
o0
,i
dV +
_
A
q
i
n
i
o0dA=0 (20)
The same procedure is applied to the balance equation (3), by multiplying by ou
i
and integrating
over the region V:
_
V
(o
i j, j
+ f
i
j u
i
)ou
i
dV =0 (21)
On the rst term of this equation we apply the divergence theorem and the straindisplacement
equation (1):
_
V
o
i j, j
ou
i
dV =
_
V
o
i j
oc
i j
dV +
_
A
o
i j
n
j
ou
i
dA (22)
Figure 1. Body V in equilibrium with mechanical and thermal mixed boundary conditions: (a) mechanical:
the body is subject to a pressure on A
o
and to an imposed displacement in A
u
. Small displacements
variations ou
i
are also represented and (b) thermal: the body is subject to a heat ux on A
q
and to a
temperature shift A
0
. Small temperature shift variations o0 are also represented.
Copyright q 2008 John Wiley & Sons, Ltd. Int. J. Numer. Meth. Engng (2008)
DOI: 10.1002/nme
A FINITE ELEMENT FORMULATION FOR THERMOELASTIC DAMPING
Then integral equation (21) may be rewritten as:

_
V
o
i j
oc
i j
dV +
_
A
o
i j
n
j
ou
i
dA+
_
V
( f
i
j u
i
)ou
i
dV =0 (23)
Subtracting Equation (23) from Equation (20) we obtain the weak form of the thermoelastic BVP:
_
V
T
0
so0dV
_
V
q
i
o0
,i
dV +
_
V
o
i j
oc
i j
dV =
_
A
q
i
n
i
o0dA
+
_
A
o
i j
n
j
ou
i
dA+
_
V
( f
i
j u
i
)ou
i
dV (24)
We note that the second term in the left-hand side of (24) is linked to the thermoelastic dissipation.
By substituting the constitutive equations (7), (8) and (9) we obtain:
T
0
_
V
[
i j
c
i j
o0dV +
_
V
c
E

0o0dV +
_
V
o0
,i
k
i j
0
, j
dV
+
_
V
oc
i j
C
i j kl
c
kl
dV
_
V
oc
i j
[
i j
0dV +
_
V
ou
i
j u
i
dV
=
_
A
q
i
n
i
o0dA+
_
A
o
i j
n
j
ou
i
dA+
_
V
f
i
ou
i
dV (25)
Appropriate boundary conditions must be specied to guarantee that the thermoelastic BVP is
well-posed and a unique solution can be obtained [18]. In Figure 1 we show the mixed mechanical
boundary conditions on A= A
o
A
u
and the mixed thermal boundary conditions on A= A

0
A
q
:
o
i j
n
j
= o on A
o
, u
i
= u on A
u
q
i
n
i
=

Q on A

Q
, 0=

0 on A

0
(26)
The thermoelastic damping is due to irreversible heat ux inside the body. In the variational
equation (25) we identify the term linked to this dissipation mechanism and consider it as the
variation of a dissipation function D:
oD=
_
V
o0
,i
k
i j
0
, j
dV (27)
The loss angle [ may be evaluated according to Equation (15).
3.2. Connections with Biots principle
Our dissipation function is closely related to the dissipation function in Biots variational principle.
According to Biot [12, 13] a variational principle can be obtained in the framework of irreversible
thermodynamics by considering (u
i
, H
i
) as state variables, where H
i
is related to the entropy of the
system and to the heat ux according to equations: s =H
i,i
, q
i
=T
0

H
i
. Consider small variations
(ou
i
, oH
i
) of the state variables in the solid V. Biot introduces two invariants, the thermoelastic
potential and the dissipation function D

, whose variation is:

o(+D

) =
_
A
0oH
i
dA+
_
A
o
i j
n
j
ou
i
dA+
_
V
( f
i
j u
i
)ou
i
dV (28)
Copyright q 2008 John Wiley & Sons, Ltd. Int. J. Numer. Meth. Engng (2008)
DOI: 10.1002/nme
E. SERRA AND M. BONALDI
The thermoelastic potential is =W+P, where Wis the isothermal strain elastic energy and P
is the heat potential. Further details can be found in [14]. The isothermal elastic energy variation
oW, the heat potential variation oPand the dissipation function variation oD

can be expressed as:

oW=
_
V
o
i j
oc
i j
dV+
_
V
0[
i j
oc
i i
dV, oP=
c
E
T
0
_
V
0o0dV, oD

=
T
0
2
_
V
k
1
i j
o(

H
j

H
i
) dV (29)
The dissipation function is related to the irreversible entropy produced inside the body due to the
volume deformations. In fact we have
_
D

dt =
_
T
0
2
_
V
k
1
i j

H
j

H
i
dVdt (30)
and taking into account the denition of

H
i
,

H
j
and the symmetry of the conductivity tensor
k
i j
=k
j i
(Onsagers reciprocity relations), we obtain on a cycle:
_
D

dt =
_
1
2T
0
_
V
0
,i
k
i j
0
, j
dVdt (31)
Comparing the denition of the dissipation function (14) it is easy to nd that:
D=D

T
0
(32)
This relation gives a deeper insight into the dissipation function D, which we use in the following
to evaluate the thermoelastic loss angle. We point out that a nite element formulation based on
Biots variational principle is also possible but it would require the use of vector variable H
i
instead of the scalar 0, increasing the size of the algebraic system of equation.
4. THERMOELASTIC ELEMENTS
In the following we dene two elements for solving problems involving either thin or thick solid
bodies. For the thin structures we develop a new thermoelastic element satisfying the Reissner
Mindlin plate theory with a through-the-thickness linear approximation of the temperature eld.
For thick bodies, we use a classical quadratic iso-parametric element.
4.1. 8-Node thermoelastic ReissnerMindlin plate element
We seek a coupled-eld element which takes into account the in-plane and out-of-plane mechanical
deformations and the thermal balance but at the same time we want to preserve a small number
of DOFs. For this reason we derive the deformation along the plate thickness due to bending
and shear stresses by following a given plate theory. A new planar 8-node serendipity coupled-
eld element is developed on the basis of the ReissnerMindlin theory. The basic hypothesis of
ReissnerMindlin theory are:
(i) the displacements in an external coordinate system must satisfy the relations:
u
1
=x
3
[
x
1
(x
1
, x
2
), u
2
=x
3
[
x
2
(x
1
, x
2
), u
3
=u
3
(x
1
, x
2
) (33)
where [
x
1
, [
x
2
are the rotations while x
3
is the element local coordinate along the plate
thickness. This implies that the normal remains straight after deformation, but it is not
Copyright q 2008 John Wiley & Sons, Ltd. Int. J. Numer. Meth. Engng (2008)
DOI: 10.1002/nme
A FINITE ELEMENT FORMULATION FOR THERMOELASTIC DAMPING
necessarily perpendicular to the deected mid-surface. The transverse displacement u
3
does
not vary along the thickness.
(ii) the stress along the plates thickness is negligible (o
33
=0).
(iii) a shear correction factor m=
5
6
is introduced to account for the parabolic distribution of
the shear stress along the thickness.
A full three-dimensional description of the thermoelastic damping requires a proper handling of the
thermal balance in the element. As shown in Equation (16), the heat ux due to the thermoelastic
coupling represents the source term in the heat equation and is related to the time derivative of
the trace of the strain tensor c
kk
=c
11
+c
22
+c
33
. According to the hypothesis (i), we can write the
in-plane component of the trace as
c
in-plane
c
11
+c
22
=x
3
_
*[
x
1
*x
1
+
*[
x
2
*x
2
_
(34)
while the out-of-plane can be derived considering hypothesis (i), (ii) and the constitutive
equation (7):
c
out-of-plane
c
33
=
z
z+2j
x
3
_
*[
x
1
*x
1
+
*[
x
2
*x
2
_
(35)
Combining the above relations and using the denitions of Lam e constants, the heat generated
by the thermoelastic coupling for small harmonic oscillations with respect to the equilibrium
temperature T
0
, may be written as
q
ThEl
=jc
E:T
0
(1v)
x
3
_
*[
x
1
*x
1
+
*[
x
2
*x
2
_
(36)
where c is the angular frequency.
Hence, the heat ux generated inside the element due to thermoelastic damping depends linearly
on the distance from the plates neutral plane and changes sign when going from the upper to the
lower surface.
Now we discuss the linear approximation of the temperature shift following the approach
developed by Lifshitz and Roukes [1], which is here extended to thin plates. We consider an
isotropic body and we analyze only temperature changes along the plate thickness. The temperature
prole can be determined solving heat equation (16) that becomes
*
2
0
*x
2
3
=
jcc
E

_
0
E:T
0
c
E
(1v)
_
*[
x
1
*x
1
+
*[
x
2
*x
2
_
x
3
_
(37)
and whose solution is
0
E:T
0
c
E
(1v)
_
*[
x
1
*x
1
+
*[
x
2
*x
2
_
x
3
= Asin(x
3
)+B cos(x
3
) (38)
Copyright q 2008 John Wiley & Sons, Ltd. Int. J. Numer. Meth. Engng (2008)
DOI: 10.1002/nme
E. SERRA AND M. BONALDI
where =(1+j)cc
E
/. Applying the boundary conditions (no heat ow across the boundaries of
the plate), we derive the constants A and B and nd the solution:
0(x
1
, x
2
, x
3
) =
E:T
0
c
E
(1v)
_
*[
x
1
*x
1
+
*[
x
2
*x
2
_
_
x
3

sin(x
3
)
cos(h/2)
_
(39)
where *[
x
1
/*x
1
+*[
x
2
/*x
2
is the sum of the middle-plane curvatures of the plate as it is shown
in Figure 3 and h is the plate thickness. Using a Taylor series expansion, Equation (39) may be
written as:
0(x
1
, x
2
, x
3
) =
E:T
0
c
E
(1v)
_
*[
x
1
*x
1
+
*[
x
2
*x
2
__
x
3
_
1
1
cos(h/2)
_
+
1
6

2
x
3
3
cos(h/2)
+o(x
3
3
)
_
(40)
If we consider only the linear part in order to maintain the equations system of small size, only
one coefcient is needed to compute the temperature shift. If we would add the cubic term, an
additional DOF in the temperature shift should also be considered. For thin plates the error made
with the rst order approximation is small as it will be shown in Section 6. Hence we can write,
the temperature shift as:
0(x
1
, x
2
, x
3
) =

0(x
1
, x
2
)
2
h
x
3
=

0(x
1
, x
2
)
3
(41)
To summarize, this element has four nodal DOFs ([
x
1
, [
x
2
, u
3
, 0), and the in-plane functions
are the well-known quadratic serendipity functions that satisfy completeness and compatibility
requirements. These functions are well suited in modeling thin or moderately thick structures as
described in Reference [19]. We assume also a linear dependence between natural coordinate
3
and the cartesian coordinate x
3
in the coordinate transformation. The element is represented in
Figure 2(a) where the transformation from the natural space to the physical space is also shown.
4.1.1. Shape functions. The rotations along the axis x
1
, x
2
and the displacement of the middle
plane u
3
inside each element are obtained by interpolating the nodal values with shape functions:
[
x
1
=
8

i =1
N
(i )
[
(i )
x
1
, [
x
2
=
8

i =1
N
(i )
[
(i )
x
2
, u
3
=
8

i =1
N
(i )
u
(i )
3
(42)
The temperature shift 0 is expressed as:
0=
3
8

i =1
N
(i )
0
(i )
(43)
The displacement vector {u} in the element and the temperature shift may be written using matrix
notation and Equations (42) and (43), respectively:
{u} = N
,
{,
e
} (44)
0 = N
0
{0
e
} (45)
Copyright q 2008 John Wiley & Sons, Ltd. Int. J. Numer. Meth. Engng (2008)
DOI: 10.1002/nme
A FINITE ELEMENT FORMULATION FOR THERMOELASTIC DAMPING
Figure 2. (a) The 8-node thermoelastic element with its node connectivity and geometry. T
1
is the transformation from the natural coordinate to cartesian coordinate and (b) the 20-node
thermoelastic element with its node connectivity and geometry. T
2
is the transformation from
the natural coordinate to cartesian coordinate.
Figure 3. Prole of the temperature shift 0 along the thickness. The dashed line is the real part of the
temperature shift along the plates thickness according to Equation (39). The continuous line is the linear
approximation used to describe the through-the-thickness temperature shift of our plate element.
where {,
e
} ={[
(1)
x
1
[
(1)
x
2
u
(1)
3
. . . [
(8)
x
1
[
(8)
x
2
u
(8)
3
}
T
and {0
e
} ={0
(1)
. . . 0
(8)
}
T
. Natural coordinates and
cartesian coordinate are related by the transformation:
x
1
=
8

i =1
N
(i )
x
(i )
1
, x
2
=
8

i =1
N
(i )
x
(i )
2
, x
3
=
h
2

3
(46)
where x
(i )
1
, x
(i )
2
are the cartesian coordinate of the i th node. Here we assume a linear dependence
of through-the-thickness coordinate x
3
and h is the plates thickness which does not involve any
nodal coordinate.
Copyright q 2008 John Wiley & Sons, Ltd. Int. J. Numer. Meth. Engng (2008)
DOI: 10.1002/nme
E. SERRA AND M. BONALDI
4.1.2. Straindisplacement and constitutive equations. Let us rst derive the straindisplacement
relation. Following the theory of plates [20], the strain tensor of a plate can be represented by ve
components and can be decoupled into bending and shear component:
{c} =
_
c
b
c
s
_
(47)
where c
b
={c
11
, c
22
, 2c
12
} are the bending components and c
s
={2c
23
, 2c
31
} are the shear compo-
nents. By using Equation (33) the bending strain tensor becomes:
{c
b
} =x
3

*[
x
1
*x
1
*[
x
2
*x
2
*[
x
1
*x
2
+
*[
x
2
*x
1

(48)
while the shear strain tensor:
{c
s
} =

*u
3
*x
2
[
x
2
*u
3
*x
1
[
x
1

(49)
From assumption (ii) the following relation between strain components can be derived:
c
33
=
z
z+2j
(c
11
+c
22
) (50)
Thus, the straindisplacement relation may be written as follows:
{c} =B
,
{,
e
} (51)
where
B
,
=
_
B
b
,
B
s
,
_
(52)
with
B
b
,
=x
3

*N
(1)
*x
1
0 0
*N
(8)
*x
1
0 0
0
*N
(1)
*x
2
0 0
*N
(8)
*x
2
0
*N
(1)
*x
2
*N
(1)
*x
1
0
*N
(8)
*x
2
*N
(8)
*x
1
0

(53)
Copyright q 2008 John Wiley & Sons, Ltd. Int. J. Numer. Meth. Engng (2008)
DOI: 10.1002/nme
A FINITE ELEMENT FORMULATION FOR THERMOELASTIC DAMPING
and
B
s
,
=

0 N
(1)
*N
(1)
*x
2
0 N
(8)
*N
(8)
*x
2
N
(1)
0
*N
(1)
*x
1
N
(8)
0
*N
(8)
*x
1

(54)
The stress constitutive equation (7) for the plate is
{o} =C
,
B
,
{,
e
}{[}
T
N
0
{0
e
} (55)
where elastic matrix can be written as:
C
,
=
_
C
b
,
0
0 C
s
,
_
=

z+2j

z 0

z

z+2j 0
0 0 j

0
0
_
j 0
0 j
_

(56)
and
{[}=

(3

z+2j):
(3

z+2j):
0
0
0

(57)
with

z=2zj/(z+2j). In the rst term of Equation (55) the stress tensor is decoupled into its
bending and shear component {o} ={o
b
o
s
} ={{o
11
o
22
o
12
}{o
23
o
31
}}. The entropy constitutive law
(8) may be written in matrix form as:
s ={[}
T
B
,
{,
e
}+
c
E
T
0
N
0
{0
e
} (58)
Fouriers law in matrix notation is
{q}=K
0
{0}=K
0
B
0
{0
e
} (59)
Copyright q 2008 John Wiley & Sons, Ltd. Int. J. Numer. Meth. Engng (2008)
DOI: 10.1002/nme
E. SERRA AND M. BONALDI
where
B
0
=J
1

3
*N
(1)
*
1

3
*N
(8)
*
1

3
*N
(1)
*
2

3
*N
(8)
*
2
N
(1)
N
(8)

(60)
and J
1
is the Jacobian matrix of the coordinate transformation.
4.2. 20-Node thermoelastic element
A three-dimensional element is needed to describe thick structures. We chose a hexahedral element
of the serendipity family containing only exterior nodes (20-node element shown in Figure 2(b)).
This element has four nodal DOFs (u
i
, 0) with i =1, . . ., 3 and a midside node along its edges.
In this case the straindisplacement and constitutive relations are formally the same as it is
described in the above section. In deriving the straindisplacement and constitutive matrices the
six component of the stress and the strain tensor must be considered. In this case the unknowns
vectors should also be changed according to the denition of the 20-node shape function as
{,
e
} ={u
(1)
1
u
(1)
2
u
(1)
3
. . . u
(20)
1
u
(20)
2
u
(20)
3
}
T
and {0
e
} ={0
(1)
. . . 0
(20)
}
T
.
5. u0-BASED FINITE ELEMENT FORMULATION
To get the numerical solution, we must rewrite the variational equation (24) in matrix notation
and in the frequency domain. From now on we consider only the case of null body forces f
i
. The
variational equation (24) becomes:
c
2
_
V
j o{u}
T
{u} dV +jc
_
V
T
0
o0sdV
_
V
o{0}
T
{q} dV +
_
V
o{c}
T
{o} dV
=
_
A
q
o0{

Q} dA+
_
A
o
o{u}
T
{ o} dA (61)
where c is the angular frequency. We point out that all unknowns are intended as phasors. The
solution in time is given by as the real part of these phasors. In order to derive the algebraic system
of equations we observe that each volume integral on the left-hand side of (61) is the sum of the
volume integral over the elements of the mesh. Each integral must be computed by substituting
the shape functions for the 8-node or 20-node thermoelastic element and the constitutive relations.
Each term of (61) may be derived separately, hence:
rst term becomes
c
2
_
V
j o{u}
T
{u} dV =c
2
{o,
g
}
T
M
,,,
{,
g
} (62)
Copyright q 2008 John Wiley & Sons, Ltd. Int. J. Numer. Meth. Engng (2008)
DOI: 10.1002/nme
A FINITE ELEMENT FORMULATION FOR THERMOELASTIC DAMPING
where
M
,,,
=

e
_
V
e
jN
T
,
N
,
dV (63)
The shape function matrix N
,
has size 3n1 while the mass matrix M
,,,
has size 3n3n where n is
the number of nodes into volume V which is discretized by the 8-node or by 20-node thermoelastic
elements. The global displacement vector is dened as: {,
g
} ={[
(1)
x
1
[
(1)
x
1
u
(1)
3
. . . [
(n)
x
1
[
(n)
x
1
u
(n)
3
}
T
or
{,
g
} ={u
(1)
1
u
(1)
2
u
(1)
3
. . . u
(n)
1
u
(n)
2
u
(n)
3
}
T
for the 8-node or the 20-node element, respectively.
The second term becomes after substituting the entropy constitutive equation (58):
jc
_
V
T
0
o0sdV =jc{o0
g
}
T
D
0,,
{,
g
}+{o0
g
}
T
D
0,0
{0
g
} (64)
where
D
0,,
=T
0

e
_
V
e
N
T
0
{[}
T
{B
,
} dV (65)
and
D
0,0
=

e
_
V
e
c
E
N
T
0
N
0
dV (66)
The damping matrix D
0,,
has size n3n and represents the effects of the strain eld in the heat
equation while the damping matrix D
0,0
has size nn. The global temperature shift is dened as
{0
g
} ={0
(1)
. . . 0
(n)
}
T
. N
0
has size n1.
The third term becomes
_
V
o{0}
T
K
0
{0} dV ={o0
g
}
T
K
0,0
{0
g
} (67)
where:
K
0,0
=

e
_
V
e
B
T
0
K
0
B
0
dV (68)
The stiffness matrix K
0,0
has size nn.
The fourth term becomes after substituting the stress constitutive equation (55):
_
V
o{c}
T
{o} dV ={o,
g
}
T
K
,,,
{o,
g
}+{o,
g
}
T
K
,,0
{0
g
} (69)
where:
K
,,,
=

e
_
V
e
B
T
,
C
,
B
,
dV (70)
and
K
,,0
=

e
_
V
e
B
T
,
{[}N
0
dV (71)
Copyright q 2008 John Wiley & Sons, Ltd. Int. J. Numer. Meth. Engng (2008)
DOI: 10.1002/nme
E. SERRA AND M. BONALDI
The stiffness matrix K
,,,
has size 3n3n. The stiffness matrix K
,,0
has size 3nn and represents
the effects of the thermal eld in the equation of motion.
Volume integrals are then evaluated in the elements natural coordinates. The matrix system
becomes:
c
2
_
M
,,,
0
0 0
__
{,
g
}
{0
g
}
_
+jc
_
0 0
D
0,,
D
0,0
__
{,
g
}
{0
g
}
_
+
_
K
,,,
K
,,0
0 K
0,0
__
{,
g
}
{0
g
}
_
=
_
{ f
,
}
{0}
_
(72)
The uniqueness of the solution must be guaranteed by a proper choice of boundary conditions on
displacements and thermal shifts as well as by the specication of thermal and mechanical loads.
The off-diagonal terms in the global damping and stiffness matrices represent the effect of the
two-way coupling due to thermoelasticity. These matrices are related according to the equation:
D
0,,
=T
0
K
T
,,0
(73)
5.1. Evaluation of the thermoelastic loss angle
We evaluate the thermoelastic loss on the basis of the dissipation function dened in Equation (14),
which in vector notation becomes:
D=
1
2

e
_
V
e
{0}
T
K
0
{0} dV (74)
In the time domain, we may write the temperature gradient vector as follows:
{0} ={|0|} cos(ct +[) (75)
which is the real part of the temperature gradient phasor obtained from the solution of the system
(72). Substituting (75) into (74) we obtain:
D=
1
2

e
_
V
e
{|0|}
T
K
0
{|0|} cos
2
(ct +[) dV (76)
According to Equation (15), to evaluate the thermoelastic damping we integrate the dissipation
function over a cycle:
_
cycle
Ddt =
_
2/c
0
Ddt =
1
4
_
2/c
0

e
_
V
e
{|0|}
T
K
0
{|0|} dVdt
+
1
4
_
2/c
0

e
_
V
e
{|0|}
T
K
0
{|0|} cos(2ct +2[) dVdt (77)
where we have used the cosine bisection formula. In Equation (77) we observe that the second
term vanishes and after some algebras it changes into:
_
cycle
Ddt =

2c

e
_
V
e
({0
e
Re
}
T
B
T
0
K
0
B
0
{0
e
Re
}+{0
e
Im
}
T
B
T
0
K
0
B
0
{0
e
Im
}) dV (78)
Copyright q 2008 John Wiley & Sons, Ltd. Int. J. Numer. Meth. Engng (2008)
DOI: 10.1002/nme
A FINITE ELEMENT FORMULATION FOR THERMOELASTIC DAMPING
The temperature shift module may be written as: |0| =
_
(0
g
Re
)
2
+(0
g
Im
)
2
. Hence, according to
Equation (15) the loss angle is given by:
[(c) =
{0
g
Re
}
T
K
0,0
{0
g
Re
}+{0
g
Im
}
T
K
0,0
{0
g
Im
}
cT
0
{,
g
Re
}
T
K
,,,
{,
g
Re
}
(79)
6. ANALYTICAL AND NUMERICAL RESULTS FOR THIN BEAMS
In this section we present the results obtained in a series of numerical tests to evaluate the
performances of the thermoelastic elements developed in this work. As test bench we consider
an isotropic thin cantilever silicon beam L =57mm long and w=10mm wide and with thickness
h =92m. Material properties of silicon are reported in Table I. The loss angle of this beam was
measured at a number of frequencies [22] and was found to be in agreement with the theoretical
values obtained by solving the two-way coupled equations for thermoelastic damping [1].
We consider two sets of boundary conditions for the beam: clampedfree and clampedclamped.
The isothermal resonant frequencies are given by the well-known solution of the eigenvalue problem
for bending vibrations of beams:
c
n
=a
2
n
h
L
2
_
E
12j
(80)
where a
n
={1.875, 4.694, 7.855, . . .} for clampedfree while a
n
={4.730, 7.853, 10.996, . . .} in the
clampedclamped case. The eigenfrequencies are summarized in Tables II and III for the rst ve
bending modes. The expected loss angle is computed by the relation [1]:
[
n
=
E:
2
T
0
C
p
_
6

2

6

3
sinh+sin
cosh+cos
_
(81)
where =h

c
n
/2,, , is the thermal diffusivity and C
p
is the volumetric specic heat at constant
pressure.
The thermoelastic damping analysis of the beam was implemented in a nite element solver
(APFEM) developed using Mathematica [23]. We performed a structuralthermal harmonic
analysis, forcing one mode at a time. Each modal shape is driven at its resonances by a proper
set of harmonic forces: the clampedfree beam is driven by a force applied on its tip while the
Table I. Material properties of Silicon (data from [21]) at 300 K.
Material properties
Youngs modulus, E 162.4 GPa
Poisson modulus, v 0.28
Density, j 2330 kgm
3
Thermal expansion coefcient, : 2.5410
6
K
1
Thermal conductivity, 145 Wm
1
K
1
Heat capacity per unit volume, c
E
/j 711 Jkg
1
K
1
Copyright q 2008 John Wiley & Sons, Ltd. Int. J. Numer. Meth. Engng (2008)
DOI: 10.1002/nme
E. SERRA AND M. BONALDI
Table II. Analytical (ANA) and numerical (FE) modal frequencies of the rst ve normal
bending modes for the thin clampedfree silicon beam.
Frequency (Hz) Frequency (Hz) Frequency (Hz)
a
n
(ANA) (FE 8-node element) (FE 20-node-element)
1.87 38.18 38.562 38.792
4.69 239.31 241.45 242.95
7.85 670.15 676.93 681.23
10.99 1313.15 1329.6 1338.3
14.14 2170.74 2203.5 2218.7
The numerical frequency is referred to planar discretization level (3, 60).
Table III. Analytical (ANA) and numerical (FE) modal frequencies of the rst ve normal
bending modes for the thin clampedclamped silicon beam.
Frequency (Hz) Frequency (Hz) Frequency (Hz)
a
n
(ANA) (FE 8-node element) (FE 20-node-element)
4.730 243.00 246.97 248.12
7.853 669.81 680.42 684.08
10.996 1313.26 1335.4 1343.4
14.137 2170.69 2211.3 2226.1
17.2788 3242.72 3309.3 3334.6
The numerical frequency is referred to planar discretization level (3, 60).
Figure 4. The beam is described by a single layer of 8-node plane elements or by four layers of 20-node
three-dimensional elements. A xed number of three elements are placed along the beams width, while
different discretization levels are used along the beams length. We show here planar discretization levels:
(3, 10), (3, 30), (3, 60). The rst index is the number of division along the beams while the second is the
number of division along the beams length.
clampedclamped beam is driven by a set of forces applied where the modal shape has its
maximum values. The solver reads geometry and boundary conditions data from the output of a
commercial grid generator, performs the global matrix assembly and solves the complex equation
system according to the formulation described in the previous sections. The loss angle is then
evaluated on the solution according to Equation (79). Meshes with different number of elements
(Figure 4) are used to test the convergence of the result. As shown in Figure 5, the convergence
is guaranteed with a number of 60 elements along the beams length.
Copyright q 2008 John Wiley & Sons, Ltd. Int. J. Numer. Meth. Engng (2008)
DOI: 10.1002/nme
A FINITE ELEMENT FORMULATION FOR THERMOELASTIC DAMPING
10 20 30 40 50 60
3.0x10
7
3.5x10
7
4.0x10
7
4.5x10
7
5.0x10
7
5.5x10
7
Convergence of the loss angle (f=38 Hz)
_LR
number of division along a plate's axis

_FE

_
F
E
Figure 5. Convergence analysis. The loss angle of the rst bending mode of the clampedfree beam
is evaluated by an FE analysis based on the 8-node plane element. The mesh is gradually rened by
increasing the number of elements along the beams length.
1E-7
1E-7
1E-6
1E-5
1E-4
38 Hz
240 Hz
671 Hz
1315 Hz
2174 Hz
_FE= _ANA

_
F
E
1E-6
1E-5
1E-4

_
F
E
_ANA (a)
Clamped-free beam Clamped- clamped beam
(b)
_ANA=_FE
3243 Hz
2171 Hz
1313 Hz
670 Hz
243 Hz
1E-6 1E-5 1E-4
_ANA
1E-6 1E-5 1E-4
Figure 6. 8-node thermoelastic element. Comparison of the loss angle computed with nite element [
FE
and the theoretical values [
ANA
for the clampedfree beam (a) and for the clampedclamped beam (b).
The continuous line represents the points where [
FE
=[
ANA
. The frequency values close to the data refer
to the modes listed in Tables II and III.
The results obtained with the 8-node thermoelastic element and the 20-node thermoelastic
element for the two different boundary conditions are summarized in Figure 6 and in Figure 7
respectively. Here the loss angles for the rst ve bending modes are compared with the theoretical
values given by Equation (81). An agreement within 5% is obtained for the majority of the modes
and for the two boundary conditions, conrming a good match with the theory. We remark that 8-
node thermoelastic plane elements give accurate results while requiring much lower computational
efforts than three-dimensional 20-node elements.
Copyright q 2008 John Wiley & Sons, Ltd. Int. J. Numer. Meth. Engng (2008)
DOI: 10.1002/nme
E. SERRA AND M. BONALDI
1E-7
1E-7
1E-6
1E-5
1E-4
1E-6
1E-6
1E-5
1E-4
(b) (a)
3335 Hz
2226 Hz
1343 Hz
684 Hz
248 Hz
Clamped-clamped beam

_
F
E
_
2219 Hz
1338Hz
681 Hz
243 Hz
39 Hz
_FE= _ANA

_
F
E
_
Clamped-free beam
_FE= _ANA
1E-5 1E-4 1E-6 1E-5 1E-4
Figure 7. 20-node thermoelastic element. Comparison of the loss angle computed with nite element [
FE
and the theoretical values [
ANA
for the clampedfree beam (a) and for the clampedclamped beam (b).
The continuous line represents the points where [
FE
=[
ANA
. The frequency values close to the data refer
to the modes listed in Tables II and III.
7. CONCLUSIONS
This paper presented a new nite elements formulation with thermoelastic capabilities. Starting
from a variational form of thermoelasticity, the thermoelastic damping is calculated from the
irreversible ow of entropy due to the thermal uxes that have originated from the volumetric strain
variations. We introduced a dissipation function D, which can be integrated over an oscillation
period to evaluate the dissipated energy. On this basis we developed a nite element framework for
the computation of the thermoelastic damping involving plates and/or massive bodies, with two
elements for solving problems involving either thin or thick solid bodies. For the thin structures
we proposed a plane ReissnerMindlin element with extended capabilities to model through-
the-thickness thermal effects. In this element out-of-plane mechanical deformation follows the
ReissnerMindlin assumptions while heat ux exchanges through-the-thickness are modeled by a
linear interpolation of the a temperature shift 0. We showed that this element allows an accurate
evaluation of the thermoelastic damping in a thin beam with small computational efforts. Thick
bodies may be also modeled by a standard quadratic iso-parametric element.
Our nite element framework will have immediate application in the evaluation of the thermoe-
lastic damping in multi-layered structures made by bonding of silicon wafers [24]. It could also
be used to model laminated composites of layered thin lms, adopted in the manufacture of some
microresonators [25] and in the production of high-reectivity mirrors [4]. Future extensions of
this work could include the treatment of a nite value of the thermal wave speed.
ACKNOWLEDGEMENTS
The authors are grateful to Ana Maria Alonso Rodriguez for the helpful discussions and suggestions
about the nite element formulation for thermoelasticity. This work was supported partially by European
Community (project ILIAS, c.n. RII3-CT-2004-506222) and by Provincia Autonoma di Trento (project
QL-Readout).
Copyright q 2008 John Wiley & Sons, Ltd. Int. J. Numer. Meth. Engng (2008)
DOI: 10.1002/nme
A FINITE ELEMENT FORMULATION FOR THERMOELASTIC DAMPING
REFERENCES
1. Lifshitz R, Roukes ML. Thermoelastic damping in micro- and nanomechanical systems. Physical Review B 2000;
68(8):56005609.
2. Norris AN. Dynamics of thermoelastic thin plates: a comparison of four theories. Journal of Thermal Stresses
2006; 29(2):169195.
3. Liu YT, Thorne KS. Thermoelastic noise and homogeneous thermal noise in nite sized gravitational-wave test
masses. Physical Review D 2000; 62(12):122002(10).
4. Fejer MM, Rowan S, Cagnoli G, Crooks DRM, Gretarsson A, Harry GM, Hough J, Penn SD, Sneddon PH,
Vyatchanin SP. Thermoelastic dissipation in inhomogeneous media: loss measurements and displacement noise in
coated test masses for interferometric gravitational wave detectors. Physical Review D 2004; 70(8):082003(19).
5. Zener C. Internal friction in solids. I Theory of internal friction in reeds. Physical Review 1937; 52(3):230235.
6. Zener C. Internal friction in solids. II General theory of thermoelastic internal friction. Physical Review 1938;
53(1):9099.
7. Zener C, Otis W, Nuckolls R. Internal friction in solids. III Experimental demonstration of thermoelastic internal
friction. Physical Review 1938; 53(1):100101.
8. Wong SJ, Fox CHJ, McWilliam S. Thermoelastic damping of the in-plane vibration of thin silicon rings. Journal
of Sound and Vibration 2006; 293(12):266285.
9. Nayfeh AH, Younis MI. Modeling and simulations of thermoelastic damping in microplates. Journal of
Micromechanics and Microengineering 2004; 14(12):17111717.
10. Tian X, Shen Y, Chen C, He T. A direct nite element method study of generalized thermoelastic problems.
International Journal of Solids and Structures 2006; 43(78):20502063.
11. Yi Y-B, Matin MA. Eigenvalue solution of thermoelastic beam resonators using a nite element analysis. Journal
of Vibration and Acoustics 1996; 129(4):478483.
12. Biot MA. Thermoelasticity and irreversible thermodynamics. Journal of Applied Physics 1956; 27(3):240253.
13. Biot MA. New thermomechanical reciprocity relations with application to thermal stresses. Journal of Aerospace
Science 1959; 26(7):401408.
14. Nowacki W. Thermoelasticity (2nd edn). Pergamon Press: Oxford, 1986.
15. Landau LD, Lifshitz EM. Theory of Elasticity (3rd edn). Pergamon Press: Oxford, 1986.
16. Kinra VK, Milligan KB. A second-law analysis of thermoelastic damping. Journal of Applied Mechanics 1994;
61(1):7176.
17. Bishop JE, Kinra VK. Thermoelastic damping of a laminated beam in exure and extension. Journal of Reinforced
Plastics and Composites 1993; 12(2):210226.
18. Atanackovic TM, Guran A. Theory of Elasticity for Scientists and Engineers (5th edn), vol. 1. Springer: Berlin,
1999.
19. Ozkul TA, Ture U. The transition from thin plates to moderately thick plates by using nite element analysis
and the shear locking problem. Thin-Walled Structures 2004; 42(10):14051430.
20. Timoshenko S, Woinowsky-Krieger S. Theory of Plates and Shells (2nd edn). McGraw-Hill: New York, 1987.
21. Material Properties Database, http://www.jahm.com [21 April 2008].
22. Reid S, Cagnoli G, Crooks DRM, Hough J, Murraya P, Rowan S, Fejer MM, Route R, Zappe S. Mechanical
dissipation in silicon exures. Physics Letters A 2006; 351(45):205311.
23. Mathematica, http://www.wolfram.com/ [21 April 2008].
24. Zendri JP, Bignotto M, Bonaldi M, Cerdonio M, Conti L, Ferrario L, Liguori N, Maraner A, Serra E, Taffarello L.
Loss budget of a setup for measuring mechanical dissipations of silicon wafers between 300 K and 4 K. Review
of Scientic Instruments 2008; 79(3):033901.1033901.12.
25. Prabhakar S, Vengallatore S. Thermoelastic damping in bilayered micromechanical beam resonators. Journal of
Micromechanics and Microengineering 2007; 17(3):532538.
Copyright q 2008 John Wiley & Sons, Ltd. Int. J. Numer. Meth. Engng (2008)
DOI: 10.1002/nme