FORMULATION OF COUPLED PROBLEMS OF

THERMOELASTICITY BY FINITE ELKMENTS

M. BALLA

Department of Technical Mechanics

Technical University, H-1521 Budapest

Received November 10, 1987

Presented by Prof. Dr. Gy. Beda

Abstract

When designing mechanical equipment, in numerous cases thermal effects have to be taken
into consideration besides mechanical ones, especially if they are transient effects. In engineering
practice these thermomechanical problems can be described by the different theories of thermo-
elasticity \vith adequate accuracy. This paper, after a brief survey of the basic equations of thermo-
elasticity, demonstrates its formulation by the finite element method.

Symbols

(Jij - linear stress tensor

Gij - linear strain tensor
u; - displacement vector
(l - mass density
F; - external force per unit mass
R - strength of the internal heat source per unit mass
qi heat flux per unit area
s - entropy per unit mass
T - absolute temperature
To - initial temperature
e - temperature difference
C;jkl - elasticity tensor
f3ij - thermoelasticity tensor
kij - thermal conductivity tensor
c - specific heat per unit mass at constant strain
I., It - Lame elastic constants
f3 - thermal modulus
Cl. - coefficient of linear thermal expansion
k - thermal conductivity
E - Young's modulus
v - Poisson's ratio
r - relaxation time
- time
r - three-dimensional space vector
71J - outward normal unit vector of the surface
qll - heat flux in normal direction
a - coefficient of convective heat transfer
ell' - surface temperature difference
B= - reference temperature difference
(:) - partial differentation with respect to time
( . ), i-partial differentation with respect to descartes coordinate Xi

1. Introduction

Mechanical equipments during operation are affected by various interactions,

the most significant being the mechanical and thermal effects. Mechanical and thermal
loads usually occur simultaneously and as a result, the displacement and temperature
fields are created in close connection with each other. The two fields have to be
defined simultaneously taking the relationship behveen them into account which
in practice proves to be rather difficult.
Thermomechanical processes are described by the basic equations of continuum-
mechanics and thermodynamics. In the solution of a variety of problems the appli-
cation of thermoelasticity proves to be efficient. The foundations or thermoelasticity
were laid by Duhamel and Neuman in the first half of the last century, widespread
interest in this field, however, has not developed till the middle of the 20th century.
The reason for the sudden growth of interest is that at that time the need for designing
equipment that can operate at very high temperatures arose almost simultaneously
in several dynamically developing areas of industry. Such areas were among others:
production of high-speed aeroplanes, design of space vehicles, rocket and jet en-
gines, technology of large turbines and the design of nuclear reactors.
Biot [1], Boley and Weiner [2], Parkus [3], Nowinski [4] and numerous other
scientists have dealt with the solution of the problems, and as a result of their work
the theory of classical linear thermoelasticity was created based on the solid founda-
tion of reversible thermodynamics. In recent years modified, generalized versions
 COUPLED PROBLD,fS OF THERMOELASTICITY 61

of the classical theory have been published (Lord and Shulman [5], Green and Lindsay
[6], Szekeres [10]).
The first part of this paper contains a brief summary of the basic equations of
classical and generalized linear thermoelasticity. The second part introduces a finite
element method which is suitable for solving two-dimensional problems.

2. Basic equations of linear thermoelasticity

2.1. Classical linear thermoelasticiry

Classical linear thermoelasticity is based on the following fundamental equa-

tions [4]:
- kinematic relation:
1
8ij = 2 (Ui.j+Uj,i)' (2.1)
- equations of motion:
uij,j+QEi = QUi' (2.2)

uij = Uji' (2.3)

- energy-scale equation:
qi,i+ Q(Tos-R) = 0, (2.4)
- constitutive equations:
Uij = CjjkZ8kZ + !3ij8, (2.5)

qj = -kij8,j, (2.6)
QC
QS = To 8-j3ij8ij. (2.7)

From the equation of motion (2.2) and the energy-scale equation (2.4) using
the linear constitutive equations (2.5)-(2.7) we get the most general basic equations
of linear thermoelasticity:
(CijkZ8kl),j+(!3ij8),j + QEi- QUi = 0, (2.8)
(k 1).. 8 ,J,
.) l.-nc(}+oR+T,oj3
~ -
.. e.. = O.
lJ I)
(2,9)
In the case of homogeneous, isotropic material with respect to (2.1), (2.8) and (2.9)
can be transcribed as:
IlUi,jj+().+P)Uk,ki+j38,i+QEi-QUi = 0, (2.10)
k8,u-Qc(}+QR+!3Toll i,i = O. (2.11)

These equations are the coupled field equations referring to Uj and 8 field va-
riables, for homogeneous and isotropic materials of the classical theory of linear
62 M. BALLA

thermoelasticity. The relationship between the displacement and temperature fields

is shown by the third term of the equation of motion and the fourth term of the
heat-conduction equation. In the absence of these the equations simplify to the
Lame equation of classical elasticity and the Fourier's heat-conduction equation
of classical thermodynamics.

2.2. Generalized linear thermoelasticity theory

From the field equations of the classical theory the equation of motion (2.10)
is hyperbolic type in case of known temperature field, so it provides finite speed for
the spread of elastic disturbances. The heat-conduction equation (2.11) based on
the classical Fourier's heat-conduction law is, however, parabolic, so the classical
theory bears on itself the paradox of infinite speed of thermal signals. To resolve
this, in the last two decades so-called generalized thermoelasticity theories were
created (e.g. Lord and Shulman [5], Green and Lindsay [6), Szekeres [10]) which
are based on hyperbolic type heat-conduction equations and these result in finite
speed of thermal signals. In these the different authors modified the heat-conduction
equation of thermoelasticity on the basis of different principles.
In the foIIowing we are going to deal with the generalized theory based on the
modification of Fourier's classical heat-conduction law. Cattaneo's and Vernotte's
modified version of the heat-conduction law is the following for anisotropic material:

( 1+'r~)
o tq.' =-k··f)
'J
'.
,j
(2.12)

If instead of equation (2.6) we regard the above equation as the constitutive equation
describing heat flux, then the fonowing equation can be derived for isotropic material:

(2.13)

This equation and the equation of motion (2.10) constitute the whole system of
field equations of generalized thermoelasticity.
The relaxation time appearing in the equation is a new material property, the
value of which according to literature is between 10- 14 and 10- 10 sec. The genera-
lized heat-conduction equation (2.13) in '0=0 case is naturaIIy simplified to the
classical heat-conduction equation (2.11).
 COUPLED PROBLEMS OF THERMOELASTICITY 63

2.3. Initial and boundary conditions

When formulating a specific practical problem, initial and boundary conditions

have to be attached to the field equations of the classical or generalized theory.
Initial conditions:
u;(r, t = 0) = uiO(r), lii(r, t = 0) = liiO(r),
8(r, t = 0) = 8o(r), Ber, t = 0) = Bo(r). (2.14)
e
The inital condition referring to has to be given only when the generalized theory
is applied.
Boundary conditions:
Let us mark the examined range with 'V', its boundary with 'A'. The boundary
conditions can be divided into two large groups: mechanical and thermal boun-
dary conditions:
1. mechanical boundary conditions:
a) kinematic boundary condition (prescribed displacement):
(2.15)
b) dynamic boundary condition:
(2.16)
2. thermal boundary conditions:
a) prescribed temperature condition:
81 A 1 = &, (2.17)
b) heat flow and convection boundary condition:
Qil1ilA, = ii,,+a(8 w -8=). (2.18)

3. Finite element analysis of two dimensional thermoelastic problems

The field equations (2.10) and (2.13) of the generalized theory of linear thermo-
elasticity make a coupled partial differential equation system consisting of four
equations. The closed form solution of this system belonging to given initial and
boundary conditions \vith the exception of a few trivial cases cannot be produced
for most practical problems. Modern computers and numerical methods have not
been at the disposal of design-engineers' until very recently. For this reason to
alleviate the problems various simplified forms of the basic equations were created
by approximate assumptions based on practical experience. When designing heavy-
duty equipment often the results gained from simplified equations are not satisfac-
64 M. BALLA

tory. In these cases the calculations have to be made on the basis of all the field
equations of thermoelasticity. This problem can usually be solved with the appli-
cation of numerical methods only. Among the different numerical methods the
finite element is probably the only one, that can most effectively be used for solving
the problems of thermoelasticity. In the following the main steps of deduction of a
finite element method suitable for solving two-dimensional problems will be outlined.
When deducing finite element schemes suitable for solving mechanical problems
the starting point is usually a variational principle. A variational principle exists
for the problems of thermoelasticity as well, which was elaborated by Keramidas and
Ting [8] through the utilization of Biot's [7] results. This principle contains the so-
called heat displacement instead of temperature which makes taking heat flow and
convection boundary conditions into account more difficult when solving practical
problems. We would like to deduce the finite element scheme referring to the field
variables of displacement and temperature, but there is no variational principle for
these variables, so the scheme \vill be produced through a special weighted residual
method, the Galerkin-method.

3.1. Discretization in space

In the case of two-dimensional problems only three scalar field functions, the
'x' directional 'u' and 'y' directional 'v' components of the displacement vector and
the temperature difference have to be defined.
For finite element formulation of the problem the so-called finite element scheme
referring to one element is needed. To deduce this the approximation of the field
functions within one element wiII be given - in accordance with the basic principles
of the method - with the local unknowns (u i , Vi' a;) and the suitably chosen shape
functions (Nt, N!):
n
uE(x, y, t) = Z UiU) Nt(x, y),
i=l

n
VE(X, y, t) = Z v i (t)N1(x, y), (3.1)
i=l

n
aE(x, y, t) = Z ai(t) Nfex, y).
i=l

Let us now assume that the nodal values are the continuous functions of time, thus
the discretization in time and space is separated. The Nt and NiB symbols signify
that various shape functions can be applied for the interpolation of displacement
and temperature fields. The n index that appears in the equations is the number of
 COUPLED PROBLEMS OF THERMOELASTICITY 65

nodal points in the applied finite element. Let us arrange the nodal unknowns of the
element into one vector:

.E ~ [::j' where dE ~ [::1' E

9 = [~ll'
(3.2)

ell
I~nl
LVII]

and the shape functions into one matrix:

N= [N 0], d
where Nd = [Ni 0 N# 0 ... N: 0] , N° = [Nf ... N%J.
o N° o Ni 0 Ng... 0 Nj~
(3.3)
The approximation of the field functions within one element with the above symbols:
(3.4)

Let us first do the discretization of the equation of motion. In order to take the
dynamical boundary condition simply into account we start out from the following
form of the equation:
(3 5)

If the approximation of the field functions (3.1) is inserted into equation (3.5), then
we get a so-called residual that usually differs from zero (from here on the invar able
symbols will be used):

[1111]
1172 a . v + QF - QU"E) ;::" 0 •
= (E" (3.6)

The weighted residual method prescribes that the weighted integral value of this
residual should be zero:
.r S·mdQ
nE
= o. (3.7)

(Q E is the domain of the finite element). If we use the shape functions as S weighted
functions, then we receive Galerkin's method:

f Nd
nE
T
• (a E • V +QF- QUE) dQ = o. (3.8)

We transform the first member that appears in the integral:

5 P. P. M. 33/1-2
 M. BALU

We transform the area integral of the first term of the equation (3:9) into a line
integral with the help of the Gauss-Ostrogradsky theorem which enables us to
take the dynamical boundary condition (2.16) into consideration:

J (Nd T
·'ITE). vdQ = J Nd T
• (lE. ndS = JNd T
• pdSp • (3.10)
nE SE s;

By mtegra:ting term by term in equation (3.8) and using equations (3.1) and (3.10)
'we get:

f (Nd T
oV):(lEdQ- f Nd T
• pdSp - f Nd T
eFdQ + f Nd T
QUE dQ = O. (3.11)
nE s; nE nE

The constitutive equation referring to the stress tensor and the approxlination of the
field functions (3.1) together with equation (3.11) give the folIowing matrix-differen-
tial equation:
(3.12)

This matrix equation is the form of equation of motion discretized in space, which
refers to one element.
The quantities appearing in the equation (3.12) are as follows:

MP.d
l] =
!"lE
f 0 [NdNt!
.... 1)

0 NfN}
0] dQ
.. , (3.l3/a)

(3. 13/b)

K6° = i j3 [a~~Lj dQ, (3. 13/c)

oNi
<:I
N~J . . = 1, ... ,11
l,j
uy

(3. 13/d)

where D1 , D 2 , D3 material properties in the case of plane strain are:

D _ E(1-v) Ev E
(3. 14/a)
1 - 0+v)(1-2v)' D2 = (1+v)(1-2v)' D3 = 2(1+v)'
 COUPLED PROBLEMS OF THERMOELASTICITY 67

in case of plane stress they are:

E Ev E
Dl = 1-v2 ' D2 = I-v2 ' Da = 2(1+v) . (3. 14/b)

Let us now consider the following form of the heat-conduction equation:

.. +ncO-nR-/3~oU"
q'J' ~ ~ J,J
= o. (3.15)

To discretize equation (3.15) similarly to the discretization of the equation of motion,

Galerkin's method is used:
JNO T . (qE. V+/2cOE-(!R-/3ToUE. V)dQ = O. (3.16)
nE

The first term can be transformed here as well:

(3.17)
The heat flow and convection boundary conditions are taken into account similarly
to the dynamical boundary condition:

J (NOT oqF). V dQ = J NBT (qE. n)dS2 =

nE Sf

= J (NOT qn dS 2 + fNBT a«()E_()",,,) dS 2 • (3.18)

Sf Sf

By doing the integration in equation (3.16) term by term and taking equations
(3.17) and (3.18), and modified Fourier's law into account we can derive the follow-
ing equation:

JNOToV)kV()EdQ+(I+-r
nE
aa) [fNO T(!cOEdQ-
tnE
JNBT(!RdQ-
QE

- fNOT/3ToUE.VdQ+ JN8 Tqn dS2+ fNOTa«()E_()",,)dS2] =0. (3.19)

nE Sf Sf

Using the approximation (3.1) from equation (3.19) we get the form of the
modified heat-conduction equation that is discretized in space:

(3.20)

M~f.'J = f-r/3T,
nE 0
[N~• aN1IN~
ax • aNJ]
ay dQ , (3.21/a)

Mj~9 = J -r(!cNfN1dQ, (3.21jb)

nE

5*
68 M. BALLA

C'Jqg = [N~
J f3T.0 , aNJ\Nq aNJ]
~. . dQ-,
,~ (3.2I/c)
nE uX u}

i,j= 1, ... ,11

CfJ= fQcNrN~dQ+ J-r:aMN~dS2' (3.21/d)

QE sf

Ki~O = J k (aN? aN~ + (iN? aNJ)dQ+ f aNrN~dS2' (3.2I/e)

nE ax ax ay ay Sf

if>? = f a8<x>NfdS2+ f Q(R+TR)N?dQ- f (c]1l+ T4n)N?dS2 • (3.21/f)

sf nE Sf
On the basis of equations (3.12) and (3.20) the matrix form referring to one
element of the coupled field equations discretized in space:

(3.22)

From the equations referring to the elements the so-called global matrix diffe-
rential equation-system of the whole structure can be constituted:

(3.23)

o is the vector of nodal unknowns,

M the generalized mass matrix, C the generalized
capacitance matrix, K the generalized stiffness/conductivity matrix and <I> the
generalized nodal force vector.

4.2. Discretization ill time

Discretization in time can be best done by three node quadratic time finite
elements. The time dependence of the nodal unknowns of division according to
space can be approximated with the help of the nodal unknowns (Ok-I, Ok' 0HI)
of the 2L1t length time finite element and the quadratic time shape functions (Mk - I ,
M k , M HI):
(3.24)

If this approximation is inserted into equation (3.23) and the weighted residual
method is applied by choosing various weight functions, then we get various time-
schemes for the discretization of equation (3.23). From these the so-called quadratic
Galerkin-scheme received through the application of the Mk+1 weighting function
 COUPLED PROBLEMS OF THERMOELASTICITY 69

shows the most favourable stability and accuracy characteristics [9]:

+ [-M - ~ C1t+ ~ K(Llt)2] 5 k- 1- ~ (Llt)2<P k_1 + ~ (Llt)2<Pk+ : (Llt)2Wk+l'

(3.25)

where <Pk - 1 , <P k , and <Pk+1 are the values of the load vector in the nodal points of the
time element. Using scheme (3.25) the solution of the matrix differential equation-
system (3.23) can be reduced to the solution of the following linear system of equa-
tions
(3.26)

The kinematic and prescribed temperature conditions have to be taken into account
when solving the so-called effective equation-system (3.26).
With the help of the quadratic Galerkin-scheme we can determine the 5k+l
value in the (k + 1) 'th time step if the value of the nodal unknowns of the (k -1)
and /.c'th steps are known. When starting out, the scheme, in case k=1, to calculate
5 2 ,5 0 and 51 are needed. 51 can be produced by the application of the various starting
methods fro111 the vectors 50 and &0 which are defined by the initial conditions. One
of these methods is the following, so-called Crank-Nicholson method:

(3.27)

From the second step on with the help of scheme (3.25) the approxinlate solution
of the matrix differential equation-system can be produced step by step.

References

1. BrOT, M. A.: Thermoelasticity and Irreversible Thermodynamics. J. App!. Phys., 27, 240-253,
(1956).
2. BOLEY, B. A. and WEINER, J. H.: Theory of Thermal Stresses. J. Wiley and Sons, New York,
(1960).
3. PARKUS, H.: Thermoelasticity. Blaisdell Pub!. Co., Waltham Massachusetts, (1968).
4. NOWINSKI, J. L.: Theory of Thermoelasticity with Applications. Sijthoff and Noordhoff Int.
Pub!', Alphen Aan Den Rijn, (1978).
5. LORD, H. W. and SHULMAN, Y.: A Generalized Dynamical Theory of Thermoelasticity. J. Meeh.
Phys. Solids, 15, 299-309, (1967).
70 M. BALLA

6. GREEN, A. E. and LlNDSAY, K. A.: Thermoelasticity. J. of Elasticity, 2, 1-7, (1972).

7. BlOT, M. A.: Variational Principles in Heat Transfer. Clarendon Press, Oxford, (1970).
8. KERAMIDAS, G. A. and TING, E. C.: A Finite Element Formulation for Thermal Stress Analysis.
Part I.: Variational Formulation. Nucl. Eng. Design, 39, 267-275, (1976).
9. STASA, F. L.: Applied Finite Element Analysis for Engineers. CBS Publishing, New York, (1985).
10. SZEKERES, A.: Equation System of Thermoelasticity Using the Modified Law of Thermal Con-
ductivity. Periodica Polytechnica, 24, 253-261, (1980).

Miha.ly BALLA H·1521, Budapest