UCD School of Electrical, Electronic and

Mechanical Engineering

Lecture notes for the course

Computational Continuum Mechanics (CCM)

Prepared by

Ismet Demirdzic
University of Sarajevo, Bosnia and Herzegovina

Alojz Ivankovic
University College Dublin, Ireland

Noel O’Dowd
University of Limerick, Ireland
 Resume
Continuum is a mathematical idealisation (representation) of real materials such as solids,
fluids and gases. Continuum mechanics is concerned with the behaviour of such continua. It
is based on laws of conservation (of mass, energy, momentum, moment of momentum and
entropy), and is aided by a mathematical apparatus specifically developed for it. With the
advent of computer technology and accompanying computational methods, Computational
Continuum Mechanics (CCM) has become an extremely powerful tool for analysing
engineering problems.

A number of numerical methods have been developed for analysing continuum mechanics
problems. These utilise conventional numerical methods such as Finite Element (FE), Finite
Volume (FV), Boundary Integral, Finite Difference, …. The FE method has been the most
widely used for the numerical analysis of solid mechanics problems, while the FV method is
the ‘number one’ method in computational fluid mechanics. The current trend is that both
methods are being developed not only for one particular discipline of continuum mechanics
but also for interdisciplinary applications.

The CCM course is designed to introduce basics of continuum mechanics, and FV and FE
methods. By doing so, it is aimed to demonstrate how fundamental laws of physics can be
turned into powerful numerical tools.

The CCM course is also planned to provide theoretical basis for practical applications of FV
and FE methods.

i
Contents

1 Introduction 1

2 Basic concepts and definitions 4

2.1 Introduction 4
2.2 Concept of continuum 4
2.3 Continuity, homogeneity and isotropy 6
2.4 Mathematical basis 6
2.5 Elements of matrix algebra 20

3 Stresses 27
3.1 Body and surface forces 27
3.2 Stress tensor 28
3.3 Principal stresses. Stress tensor’s invariants 31
3.3 Spherical and deviatoric stresses 32

4 Deformation and flow 34

4.1 Introduction 34
4.2 Deformation 37
4.3 Motion and flow 45

5 Fundamental laws of continuum mechanics 52

5.1 Introduction 52
5.2 Fundamental laws of continuum mechanics in various forms 54
5.3 Equations for large deformations 65
5.4 Variational methods 69

6 Constitutive relationships 73
6.1 Introduction 73
6.2 Classical constitutive relations and equations of state 75
6.3 Elastic solids 78
6.4 Ideal and Newtonian fluids 81
6.5 Fourier’s law and the caloric equation of state 83

ii
7 Mathematical models 85
7.1 Introduction 85
7.2 Linear elastic solids 85
7.3 Newtonian fluids 88
7.4 Initial and boundary conditions 90
7.5 Summary 95

8 Finite Volume discretisation 96

8.1 Main Concepts and Ideas of the Finite Volume Method 98
8.2 Two-dimensional FV discretisation 124
8.3 Solution procedure 145
8.4 Summary 149
8.5 Three-dimensional problems in domains of arbitrary shape 150

9 An introduction to Finite Element method 157

9.1 Background 157
9.2 Description of boundary value problem 158
9.3 Principle of virtual work 158
9.4 Finite Element discretisation 160
9.5 Linear elastic finite element model 161
9.6 The element and global stiffness matrices 163
9.7 Numerical quadrature 170
9.8 Evaluation of stiffness matrix 173
9.9 Mapping of elements to and from the parent space 173
9.10 Computation of the nodal force matrix 178
9.11 Assembling the global stiffness and force matrices 180
9.12 Application of displacement boundary conditions and solution of the finite
element equations 181
9.13 Steps in a static linear finite element analysis 181

iii
 1. INTRODUCTION

Chapter 1

Introduction
It is important to recall that much of the technical knowledge results from only a few basic laws
of physics. They form the basis for specialised engineering disciplines such as solid mechanics,
fluid mechanics, heat transfer, etc.. Continuum Mechanics deals with these principles which are
common for practical engineering materials.

Continuum is a mathematical idealisation (representation) of real materials such as solids, fluids

and gases. Continuum mechanics is concerned with the behaviour of such continua. It is based on
laws of conservation (of mass, energy, momentum, moment of momentum and entropy), and is
aided by a mathematical apparatus specifically developed for it. While these laws are universal,
additional or constitutive laws are required to define specific materials. With the advent of
computer technology and accompanying computational methods, it has become an extremely
powerful tool for analysing vast variety of engineering problems.

A number of numerical methods have been developed for analysing continuum mechanics
problems. These utilise conventional numerical methods such as Finite Element (FE), Finite
Volume (FV), Boundary Integral, Finite Difference, etc. The FE method has been the most
widely used for the numerical analysis of solid mechanics problems, while the FV method is well
known in computational fluid mechanics. Computational Continuum Mechanics (CCM) course is
designed to introduce basics of continuum mechanics. Is not intended to be a review of numerical
methods but to briefly introduce basic principles of FV and FE methods. By doing so, it is aimed
to demonstrate how fundamental laws of physics can be turned into powerful numerical tools.

Solution of complex real-life engineering problems has become an everyday reality. Demands in
this field are on ever-increasing path. Although quite powerful, the present day numerical
methods are far from being perfect, for a number of reasons, and much effort is being put into
improving the existing methods and developing new ones. The current trend is that both FV and
FE methods are being developed not only for one particular discipline of continuum mechanics
but for interdisciplinary (‘multiphysics’) applications. The direct reliance of the FV method on
the basic physical principles, and its inherent iterativeness makes it particularly suitable for these
applications.
1
 1. INTRODUCTION
There are several steps that have to be taken, and a number of important decisions to be made
along the way, in order to solve a practical engineering problem (Fig. 1.1). Each of them
inevitably involves a number of assumptions and approximations, which all together influence
the final solution of the problem. We shall now try to briefly address the essential ones.
1) The first step is a definition of the problem at hand. Its success is often limited by our ability
to grasp the full complexity of the problem, to identify the type of materials involved, to
accurately specify initial and boundary conditions etc.

2) As the next step, one has to choose a mathematical model, which can simulate the behaviour
of the problem. There, we can distinguish between:
2.1) Problems for which an adequate mathematical description can be written (e.g. heat
conduction, elastic stress analysis, laminar fluid flow), and
2.2) Problems for which an adequate mathematical description has not yet been found (non-
linear stress analysis, turbulent fluid flow).
Of course, the group into which a given problem falls will depend on what we are prepared
to consider as an 'adequate' description.
world
(1) Real

Solution domain and ENGINEERING Material

initial and boundary PROBLEM properties
conditions
(2) Physics

Constitutive MATHEMATICAL Governing

relationships MODEL equations

SYSTEM OF
Time and space ALGEBRAIC Equations
discretisation
(3) Numerical

discretisation (DISCRETISED)
method

EQUATIONS

Linearisation of SOLUTION Linear equation

equations ALGORITHM solver
software
analysis
(4) Stress

Pre- and post- COMPUTER Numerical

processors PROGRAM algorithm

SOLUTION

Fig. 1.1 Basic steps of numerical simulation of an engineering problem

2
 1. INTRODUCTION

3) The mathematical description of continuum mechanics problems is very rarely amenable to a

closed-form analytical solution and a numerical procedure is the only alternative (except for
direct measurements of quantities of interest on a physical model or on the actual object,
which has its own problems). All numerical methods consist of transforming the
mathematical model into a system of algebraic equations. In order to achieve this, the time,
space and equations have to be discretised. In doing that a number of approximations are
made: continuum is replaced by a finite set of computational points placed in time and space,
continuous functions representing the exact solution to the mathematical model are
approximated by polynomials of a finite order (much lower than the number of
computational points, except in case of spectral methods) etc.

4) The resulting set of algebraic equations is then solved by (approximate) iterative methods, by
necessity if the equations are non-linear or the number of equations is very large and the
direct solution exceeds available computer resources, or by choice if the system is linear and
not too large. Here, we should not forget the finite arithmetics of the computers and
accompanying round-off error propagation which makes theoretically exact methods at best
approximate and at worst useless.

5) Finally, the numerical method has to be translated into a computer language and computer
program assembled. In doing this a great deal of care has to be taken regarding the suitable
data structure, the program modularity etc. Also, pre- and post- processing facilities are
important components governing utility of stress analysis software.

An introduction to the continuum mechanics is presented in Chapters 2 to 7 including concept of

continuum, field theory, forces acting on and resulting motion of a continuum control volume,
conservation laws, basic constitutive laws and mathematical models. Chapters 8 and 9 deal with
principles of FV and FE discretisations, respectively.

3
 2. BASIC CONCEPTS AND DEFINITIONS

Chapter 2

Basic concepts and definitions

2.1 Introduction
There are in essence two methods for analysing material’s motion. The first one treats a material
as a set of molecules, and the macroscopic phenomena are explained as a consequence of
molecular activity, and calculated by employing the laws of mechanics and probability (the
statistical approach). The alternative method uses the concept of continuum, which neglects a
discrete structure of matter and assumes that a material fills the space continuously. The concept
of continuum, in conjunction with phenomenological macroscopic theories based on general
laws confirmed by experiments, provide effective means for the solution of practical problems
(the phenomenological approach).

2.2 Concept of continuum

Mechanics is concerned with forces acting on a body and its resulting motion and deformation.
It is based on the concepts of time, space, matter, force and energy. It can be divided into two
broad categories:

1) Classical mechanics – it includes particles and rigid body mechanics.

2) Continuum mechanics – it includes solid and fluid mechanics.

The word continuum comes from mathematics; e.g. the system of real numbers is continuous
since there is always a real number between any two different real numbers (this means that the
system of real numbers is infinite). Time can be described by one system of real numbers, while
for 3D space three systems are required.
The concept of continuum applied to a matter is best illustrated on considering the mass density.
Let’s assume that a material occupies the volume V0 (Fig. 2.1). Let’s also consider a point P in
V0 and in a series of progressively smaller sub-volumes
V1, V2 , ...(Vn 1  Vn , P Vn , n = 1,2,...) all containing point P. If volume Vn contains mass
M n , then the mass density in P is defined as:

4
 2. BASIC CONCEPTS AND DEFINITIONS

Mn
(P) = lim . (2.1)
n Vn
Vn 0
If the mass density is defined in every point of V0 and in every time instant, then the mass is
said to be continually distributed, and
 = (x, y,z, t) , (2.2)
is a single-valued function of space and time. Here, x, y and z are Cartesian co-ordinates and t is
time. In a similar way one can define energy density, momentum density etc.

V2

V1
V0

Figure 2.1: Subsets of volume V0 containing point P

However, the real substances are made of elementary particles (molecules) and the continuum
idealisation is only applicable under certain conditions. If, for example, one attempts to define
the mass density of water the problems will occur when the size of Vn becomes comparable to
the molecular dimension of water. When Vn1 and Vn differ by one molecular dimension, then
the difference between Mn 1 / Vn 1 and Mn / Vn can be finite (or even infinite). As the particles
moves in space, the ratio Mn / Vn either does not exist or it fluctuates in time and space (Fig.
2.2).

M

3 V

Figure 2.2: Definition of continuum mass density

5
 2. BASIC CONCEPTS AND DEFINITIONS

To solve the problem, let’s allow Vn to decrease but to always remain big enough so it contains
a large number of particles. If under this condition the ratio Mn / Vn tends to a final value (P),
then this value is taken as the mass density at point P. It is usually taken that:
Mn
(P) = lim 3 , (2.3)
Vn  Vn
where  is of the order of the mean free molecular path. If this value is extrapolated to Vn  0,
then we have defined ‘mathematical’ continuum with mass density (P). In practice, this
extrapolation is usually never required. For example, the molecular dimension of water is
approximately 10-10 m, and water can be considered as a continuum in all flow problems where
any characteristic size is greater than 10-10 m. The mean free molecular path in air is ~ 5 10-8 m.

As shown above, the concept of continuum requires that the mean free molecular path (a mean
distance travelled by a molecule in between two collisions) is small in comparison with the
smallest characteristic length of the body. It means it can be applied to problems where the
particular nature of the material microstructure can be neglected. This condition is, however,
fulfilled for the majority of practical materials. The matter on which the concept of continuum
can be applied is called continuum. The part of physics dealing with motion of continuum is
called continuum mechanics.

2.3 Continuity, homogeneity and isotropy

Continuum mechanics uses three basic independent assumptions regarding the material
structure:
1) Continuity – material fills the space continuously without pores or holes being present.
2) Homogeneity – material properties are the same at all points. If the properties are different
at different points such material is called heterogeneous.
3) Isotropy – material property is the same in all directions. If the properties are different in
different directions such material is called anisotropic.

For example, atmospheric air is continuous but not homogeneous as its density varies with
altitude. Wood can be considered as continuous and homogeneous but is not isotropic. The first
assumption is a requirement for the use of continuum mechanics. The other two are not.

2.4 Mathematical basis

The concept of a continuum permits the properties of continuum, like density, displacements,
temperature etc., to be expressed as continuous (in a mathematical sense) single-valued
functions of space and time. The mathematical discipline which deals with such functions is
called the field theory (which owes a great deal of its development to the scientists' interest in
problems of continuum mechanics).

6
 2. BASIC CONCEPTS AND DEFINITIONS

2.4.1 Definition of field and notation

The field is defined as a continuous single-valued distribution of a certain quantity in space and
time. The field of quantity F, for example, will be denoted as

F = F(P,t) = F(x,t), (2.4)

where x is the position vector of point P, and t is the time.

Quantity F may be a scalar, a vector or, in general a tensor, and the respective fields are called
scalar, vector or tensor fields. A scalar is a quantity defined by 30=1 number, vector is defined
by 31=3 numbers, while quantity requiring 3n numbers is called a tensor, where n is tensor order
(also called rank). In fact, scalars and vectors can be regarded as tensors of the zero and first
order, respectively. In what follows, we shall mainly consider tensors of the second order
(defined by 9 numbers), and unless stated otherwise, the word tensor will be used to denote a
tensor of the second order.

2.4.2 Invariance of geometrical and physical quantities

Geometrical and physical quantities are normally analysed using an appropriate coordinate
system. However, their properties do not depend on the choice of the coordinate system, i.e.
they are invariant of the coordinate transformation. For example an arc length or the surface of a
cylinder are independent of the coordinate system used to describe the arc or the cylinder.
Similarly, motion does not depend on the coordinate system used to define forces acting on the
body, its path, velocity or acceleration.
Scalars, vectors and tensors exist as mathematical quantities independently of coordinate
system, but they can be defined or specified in relation to a particular system via their
components. Vector has three components while tensor has nine. The choice of the coordinate
system is arbitrary but it can considerably affect the effort in problem solving.

2.4.3 Cartesian coordinate system

One can notice that expression (2.4) is independent of the choice of coordinate system.
However, when it comes to practical calculation, one has to adopt a coordinate system and
calculate by using components of vectors and tensors in that coordinate system. The simplest
and often the most appropriate system is Cartesian orthogonal system. In what follows, the
corresponding expressions in a Cartesian coordinate system will be given alongside those in
coordinate-free (invariant) form. In such a coordinate frame the position vector is
x = x i + y j + z k, (2.5)

where x, y and z are the Cartesian coordinates and i, j and k are the unit base vectors in direction
of x, y and z coordinate axes, respectively. Thus, the expression (2.4) will also be written in the
following form

7
 2. BASIC CONCEPTS AND DEFINITIONS

F = F(x, y,z, t) , (2.6)

However, one should keep in mind that a Cartesian coordinate system, although the simplest,
may not always be the most convenient one, and some problems might be easier solved in
another coordinate system (e.g. cylindrical).

The density field

 =(x,t) or  = (x,y,z,t) (2.7)
is an example of a scalar field, the displacement field
u = u(x,t) or u = u(x,y,z,t) (2.8)
is an example of a vector field, while the stress field
 =  (x,t) or  =  (x,y,z,t) (2.9)
is an example of a tensor field.
In what follows, the scalar field
s = s (x,t) or s = s(x,y,z,t) (2.10)
the vector field
v = v(x,t) or v = v(x,y,z,t) (2.11)
and the tensor field
T = T (x,t) or T = T (x,y,z,t) (2.12)
will be considered.

Each vector field can be associated with the following three scalar fields
vx = vx (x,y,z,t),
vy = vy (x,y,z,t), (2.13)
vz = vz (x,y,z,t),
where vx , vy and vz are Cartesian components of vector v and
v = vx i + vy j +vz k. (2.14)

Similarly, each tensor field T can be associated with nine scalar fields

Txx = Txx(x,y,z,t) Txy = Txy(x,y,z,t) Txz = Txz(x,y,z,t)

Tyx = Tyx(x,y,z,t) Tyy = Tyy(x,y,z,t) Tyz = Tyz(x,y,z,t) (2.15)
Tzx = Tzx(x,y,z,t) Tzy = Tzy(x,y,z,t) Tzz = Tzz(x,y,z,t)

where Txx , Txy, ..., Tzz are Cartesian components of tensor T and

T = Txx ii + Txy ij + Txz ik +

Tyx ji + Tyy jj + Tyz jk + (2.16)

Tzx ki + Tzy kj + Tzz kk.

8
 2. BASIC CONCEPTS AND DEFINITIONS

Quantities ii, ij, ..., kk, which are sometimes written as ii, ij, ..., kk, are nine linearly
independent tensorial or dyadic products of Cartesian unit vectors and have the following
properties
ij  ji

i(a j + b k) = a ij + b ik, (2.17)

where a and b are arbitrary scalars.

2.4.4 Summation convention and index notation

Apart from writing in coordinate-free and fully expanded form, a compact index (or tensor)
notation is often used as it results in shorter expressions. Here, numerical indices are used; for
example, v1, v2, v3 instead of vx, vy, vz. Coordinates and base vectors x, y, z and i, j, k are denoted
as x1, x2, x3 and i1, i2, i3, respectively. Then, expression (2.14) can be written as:
3
v = v i i i, meaning v =  vi ii (2.18)
i =1

and expression (2.16) as:

3 3
T = Tij iiij, meaning T =   Tij ii i j . (2.19)
i =1 j =1

In doing this the standard (Einstein) summation convention is employed: whenever an index
appears twice in a monomial, it implies a summation over the range of its values (from 1 to 3 in
the present case).
Vectors and tensors are often identified by their components, and referred to as vector vi and
tensor Tij. This way of writing is called tensor or index notation. According to the convention,
an index may appear in a monomial once or twice. If an index appears once, it may take any one
of the three values 1, 2 or 3 and is called unrepeated or free index. Number of free indices
defines the quantity; 0 – scalar, 1 – vector and 2 – tensor. If an index appears twice in a
monomial, than the summation convention is applied. Such an index is referred to as a dummy
index since it must assume all three values (if replaced by another index, the meaning of the
expression remains unchanged).

2.4.5 Matrix notation

Vector and tensor components can also be written in a matrix form, which enables the matrix
calculus and terminology to be directly applied to operations with vectors and tensors. To avoid
misunderstandings, square brackets are used to denote matrices. Components of vector v can be
written in the form of a matrix with one row
v = [v] = [v x , v y , v z ] , (2.20)

9
 2. BASIC CONCEPTS AND DEFINITIONS

or with one column matrix (in the so-called transposed form)

v x 
 
vT = [v]T = v y  . (2.21)
 
v z 

Tensor components are often arranged in the form of a square matrix

Txx Txy Txz 
 
T = [T ] = Tyx Tyy Tyz  . (2.22)
 
Tzx Tzy Tzz 

2.4.6 Tensor’s definition

It was stated in Section 2.4.2 that properties of scalar, vector and tensor fields are invariant of
coordinate system transformation. In fact, the coordinate transformation rule can be used to
define the tensor.

We shall only consider the transformation of Cartesian coordinates, and so-called Cartesian
tensors. In case of arbitrary curvilinear coordinate transformation, the corresponding tensors are
called general tensors. In what follows and unless stated otherwise, the word tensor will be used
to denote a Cartesian tensor.

Coordinates’ transformation
' ' '
Let’s consider two Cartesian coordinate systems Ox1x 2 x 3 and Ox1x 2 x 3 having the same pole
O. Cosines of angles between the coordinate axes of the two systems are usually denoted as aij
'
meaning aij = cos(xi , x j ) according to:

 x1 x 2 x 3 
 '   a11 a12 a13 
 x1 a11 a12 a13   
A= '  or A =  a21 a22 a23  . (2.23)
x 2 a21 a22 a23   a31
 '  a32 a33 
 x 3 a31 a32 a33 

Position vector of an arbitrary point P can be written as

' ' ' ' ' '
x = x1i1 + x 2i 2 + x 3i 3 = x1i1 + x 2i 2 + x 3i 3 , (2.24)

and the relation between the coordinates of the point as

' ' ' T '
xi = aij x j  x i = a ji x j or x = Ax  x = A x . (2.25)

Matrix A is called rotation matrix and it can be shown it is orthogonal, i.e. A-1=AT. Also,
det(A)=1. Transformation (2.25) is called orthogonal transformation.

10
 2. BASIC CONCEPTS AND DEFINITIONS

Formal definition of tensors

Since the coordinates of an arbitrary vector can be defined by the coordinates of its starting and
ending points, then:
v 'i = aij v j  vi = a ji v'j or v' = Av  v = AT v' . (2.26)

It follows that coordinates of a dyad uv can be transformed as:

' '
ui v j = (aipu p )(a jq vq ) = aip a jq u pv q . (2.27)

The expression (2.27) can be generalised giving an equation for the transformation of tensor T
' ' T
Tij = aip a jq Tpq or T = A TA . (2.28)

The following formal mathematical definitions can be given by comparing the transformation
laws for tensors of the order 0, 1 and 2 with the coordinate transformation law (2.25):
1) A scalar is a quantity which has only one component invariant of the coordinate system
transformation.
2) A vector is a quantity defined by three components which transform in the same way as the
coordinates of the position vector (compare (2.25) and (2.26)).
3) A tensor (of the second order) is a quantity defined by nine components which transform
according to (2.28). Higher order tensors can be defined in a similar way.

Alternative definition of tensors

An alternative way of defining a second order tensor is by using it as a linear vector operator,
which relates a vector u with vector v in the following way
v = Tu or vi = Tij u j or v = Tu . (2.29)

Let’s demonstrate that these two definitions are identical. It follows from (2.26) and (2.29) that
' ' '
vi = aipv p = aipTpq uq = aipTpq a jqu j = aipa jqTpq u j , (2.30)

and from (2.29) that

' ' '
v i = Tij u j . (2.31)

By comparing expressions (2.30) and (2.31) one gets

'
Tij = aip a jq Tpq , (2.32)
which is equivalent to (2.28).

11
 2. BASIC CONCEPTS AND DEFINITIONS

2.4.7 Algebraic operations

Scalar multiplication

The product of a scalar s and a vector v is a vector

sv = sviii = (svx, svy, svz), (2.33)

and the product of a scalar s and a tensor T is a tensor

sT xx sTxy sTxz 

 
sT = sTij iiij = sT yx sTyy sTyz  . (2.34)
sT sT sT 
 zx zy zz 

Addition and subtraction

The sum of two vectors a and b is a vector

a + b = (ai + bi)ii = (ax + bx, ay + by, az + bz), (2.35)

and the sum of two tensors S and T is a tensor

S xx + Txx S xy + Txy S xz + Txz 
 
S + T = (Sij + Tij)iiij = S yx + Tyx S yy + Tyy S yz + Tyz  . (2.36)
S + T S + T S + T 
 zx zx zy zy zz zz 

Multiplication of two vectors

In relation to vectors one can define the following multiplication operations:

1) Scalar (inner, dot) product of two vectors, whose result is a scalar
a b = ai bi = ax bx + ay by + az bz. (2.37)
2) Vector (cross) product of two vectors, whose result is a vector
i j k
a x b = eijk aj bk ii = ax ay az , (2.38)
bx by bz

where eijk is a permutation (Levy-Civita) symbol defined as

 1 if ijk is an even permutation of indices (123, 231, 312);
eijk =(ii  i j ) i k = 1 if ijk is an odd permutation of indices (132, 321, 213); (2.39)
 0 if any two indices are equal to each other.

12
 2. BASIC CONCEPTS AND DEFINITIONS

3) Tensor (dyadic) product of two vectors, whose result is a tensor (dyad)

ax bx ii + ax by ij + ax bzik + a x bx a x by a x bz 
 
ab = aibj iiij = ay bx ji + a y by jj + a y bz jk + = a y bx a y by a y bz  . (2.40)
az bx ki + azby kj + az bz kk a b a b a b 
 z x z y z z 

Multiplication of a vector and a tensor

Here, the following multiplication operations can be defined:
1) Scalar (inner, dot) product of a vector v and a tensor T, whose result is a vector

(v x Txx + v y Tyx + vz Tzx )i +  v x Txx + v y Tyx + v z Tzx 

 
v T = vi Tij ij = (v x Txy + v y Tyy + vz Tzy )j + =  v x Txy + v y Tyy + v z Tzy  . (2.41)
(v x Txz + v y Tyz + vz Tzz )k v T + v T + v T 
 x xz y yz z zz 

2) Scalar product of a vector v and a dyad ab, whose result is a vector

v (ab) = vi ai b ji j , (2.42)

and the following is valid

v (ab) = (v a)b . (2.43)

3) Vector product of a vector v and a tensor T, whose result is a tensor

v  T = eijkv iT jl i k i l . (2.44)

Multiplication of two tensors

In relation to tensors the following multiplication operations can be defined:

1) Scalar (inner, double-dot) product of two tensors, whose result is a scalar
S:T = Sij Tji = SxxTxx + SxyTyx + SxzTzx +

SyxTxy + SyyTyy + SyzTzy + (2.45)

SzxTxz + SzyTyz + SzzTzz.

2) Vector product of two tensors, whose result is a vector

S  T = eijk S jl Tlki i . (2.46)

3) Scalar product of two tensors, whose result is a tensor

S T = SikTkj iiij. (2.47)

13
 2. BASIC CONCEPTS AND DEFINITIONS

Identity tensor

Tensor I for which

v = I v = v I , (2.48)

is called identity tensor and is equal to:

1 0 0 
 
 =ij iiij = 0 1 0 . (2.49)
0 0 1

Symbol ij is known as the Kronecker delta, and is defined as

1 if i = j
ij = i i i j =  . (2.50)
0 if i  j

Trace and determinant of a tensor

The sum of the diagonal elements of the tensor T, for which the symbol tr is used, is termed the
trace of the tensor T
tr T = Tii = Txx + Tyy + Tzz, (2.51)

while the symbol det denotes the determinant of the tensor T, which is defined as
Txx Txy Txz
det (T) = det (Tij) = Tyx Tyy Tyz . (2.52)
Tzx Tzy Tzz

Symmetric and anti-symmetric tensor

A tensor whose elements are symmetric (equal to each other) with respect to the main diagonal
of the matrix (2.22), i.e.
T = TT or Tij = Tji (i, j = 1,2,3) (2.53)
is called a symmetric tensor, and a tensor where
T = – TT or Tij =  Tji (i, j = 1,2,3) (2.54)

is called an anti-symmetric (or skew-symmetric) tensor. Tensor TT = Tjiiiij is called the

transpose of T. It is important to note that every tensor may be expressed as a sum of one
symmetric and one anti-symmetric tensor in one and only one way
1 T 1 T 1 1
T = (T + T )+ (T  T )= [ (Tij + Tji )+ (Tij  Tji ) ] iiij. (2.55)
2 2 2 2
s
Ts T a Tij Tija

14
 2. BASIC CONCEPTS AND DEFINITIONS

Spherical tensor and tensor deviator

It is sometimes convenient to split a tensor into two parts, one called the spherical tensor and the
other the tensor deviator. This is done in the following way
T = Ts + Td = Tm + (T – Tm ) = [Tmij + (Tij – Tmij)] iiij =
(2.56)
Tm 0 0  Txx  Tm Txy Txz 
   
 0 Tm 0  +  Tyx Tyy  Tm Tyz  ,
 0 0 Tm   Tzx Tzy Tzz  Tm 
Spherical tensor Tensor deviator
where
1 1 1
Tm = tr T = Tii = (Txx + Tyy + Tzz ) . (2.57)
3 3 3

Vector and tensor invariants

A vector has only one invariant, its length or intensity
2 2 2
|v| = v v = v iv i = v x + v y + v z . (2.58)

A (symmetric) tensor, however, has three invariants. They are related to the so-called
characteristic equation
det (T – T I) = det (Tij –  ij) = 0 (2.59)
which can be written in the (Hamilton-Cayleyev) form
3 – IT2 – IIT – IIIT = 0, (2.60)
where
IT = tr T= Tii,
1 1
IIT = (T : T  tr 2 T) = [Tij Tji  (Tii )2 ], (2.61)
2 2
1
IIIT = det (T )=  ijk  pqr TipTjq Tkr
6
are the first, second and third invariant of the tensor T.

Principal tensor components and principal directions

The roots of eq. (2.59) or (2.60) T(1), T(2), T(3) are termed the principal tensor components. For
each T(k) (k = 1, 2,3) equation

(T  T(k) I) l (k) = 0 or (Tij  T(k) ij )l (k )

j = 0 (k = 1,2,3; no summation on k) (2.62)

15
 2. BASIC CONCEPTS AND DEFINITIONS

has non-trivial roots (l1(k) ,l 2(k ) ,l3(k) ) which represent components of the three principal
directions. If coordinate axes are taken along the principal directions, the tensor T becomes a
diagonal tensor
T(1) 0 0 
 
T =  0 T(2) 0  . (2.63)
 
 0 0 T(3) 

Inverse and orthogonal tensors

Tensor S for which

S T = T S = I or Sik Tkj = Tik Skj =  ij . (2.64)

is called inverse tensor of tensor T and is denoted as S = T-1. Tensor T is called direct tensor of
tensor T-1. The sufficient condition for the existence of an inverse tensor is det(T)0. It can be
shown that the principal components of an inverse tensor are equal to the inverse of the principal
components of the direct tensor, while their principal directions coincide.

Tensor T for which T-1= TT , i.e. for which

T T
T T = T T = I or Tik Tjk = Tki Tkj =  ij . (2.65)

is called orthogonal tensor.

2.4.8 Differential operations

Gradient of a scalar field

Each scalar field can be associated with a vector field of its gradient
s s s s
grad s = i = i + j + k = s, j .
x j j x y z
( ) (2.66)

Scalar multiplication of gradient by unit vector lo = l/|l| gives derivative of the scalar field s in
direction l defined by vector l
s
= grad s lo. (2.67)
l
Divergence, rotor and gradient of a vector field

Each vector field can be associated with a scalar field of its divergence
v j v x v y vz
div v = =
x j x
+
y
+
z
( )
= v j,j , (2.68)

16
 2. BASIC CONCEPTS AND DEFINITIONS

with a vector field of its curl:

i j k
v   
curl v = ijk k ii = , (2.69)
x j x y z
vx vy vz

as well as a tensor field of its gradient

v v v v j
grad v = i+ j+ k = ii =
x y z x i i j  v x v y v z 
 
v x v y v  x x x 
ii + ij + z ik +
x x x  v v y v z 
, or grad v =  x . (2.70)
v x v y v z  y y y 
ji + jj + jk +  v
y y y v y v z 
 x 
v x v y v  z z z 
z
ki +
z
kj + z kk = v j ,i
z
( )

Divergence of a tensor field

Each tensor field can be associated with vector field of its divergence, which has an important
application in continuum mechanics
 T xx Tyx Tzx 
 + + i +  Txx T yx Tzx 
 x y z   x + y + z 
 
T ji  T xy Tyy Tzy   Txy T yy Tzy 
div T =
x j
ii = 
 x
+
y
+
z 
j+ =
x
+
y
+
z 
( )
= T ji, j . (2.71)
 
 T xz T yz Tzz   Txz + Tyz + Tzz 
 + + k  x y z 
 x y z  

Especially, divergence of a dyad is the following vector

 (ax bx ) (ay bx ) (az bx )
 + + i +  (ax bx ) (ay bx ) (az bx ) 
 x y z   x + y + z 
 
(a j bi )  (ax by ) (ay by ) (az by )  (ax by ) (ay by ) (az by ) 
i = + + j+=  + +
z 
div (ab)=
z 
.(2.72)
x j i  x y

x y

 (ax bz ) (ay bz ) (az bz )  (ax bz ) + (ay bz ) + (az bz ) 
 + + k  x y z 
 x y z  

Del (nabla) operator

The above definitions of gradient, divergence and curl can be written in a common way using
the so called Hamilton (del or nabla) operator which is defined as

17
 2. BASIC CONCEPTS AND DEFINITIONS

1
(...) = lim
V 0 V
 n(...)dS , (2.73)
S

where V is the volume bounded by a surface S and n is the outward normal to the surface S. The
quantity (…) stands for either scalar, vector or tensor and is multiplied with n in a certain way.
For the Cartesian coordinates:
(...) (...) (...) (...)
(...) = ij = i+ j+ k. (2.74)
x j x y z
This operator has the properties of both a vector and a linear differential operator. It follows
that:
grad s = s, grad v =  v,

div v =   v, div T =   T, (2.75)

curl v =  x v.

Some useful relationships

Some useful relationships involving operators grad, div and curl are given below:
grad f (s) = f ' (s) grad s = f ' (s) s
div (sv) = s div v + v grad s = s  v + v s
div (T v) = div T  v + T:grad v = (  T)  v + T:(v)
(2.76)
div (ab) = b div a + a grad b = b( a) + a (b)
curl (sv) = s curl v + grad s  v = s  v + s  v
curl (a  b) = b grad a + a div b  a grad b  b div a

Especially:
div x = 3
curl x = 0
(2.77)
grad x = I
grad x = x
where x is a position vector.

2.4.9 Integral theorems

In the field theory flux of vector field through a surface S is defined as
V̇ =  v n dS =  v dS , (2.78)
S S

where n is the unit vector of the outer normal to the surface S, and circulation of the vector field
along a closed curve C as
 =  v dx, where dx = dxi + dyj+ dzk , (2.79)
C

18
 2. BASIC CONCEPTS AND DEFINITIONS

In connection with the above quantities one can prove the following two very important
theorems:
1) Gauss' theorem
Gauss' (Divergence) theorem relates the source of a field in a region V with the flux of that
field through the boundary surface of the region S

 v n dS =  div v dV , or more generally  n dS =   dV , (2.80)

S V S V

where  = s, v or T, and  stands for s,   v or   T, respectively.

2) Stokes' theorem
According to Stokes' theorem, the circulation of a vector field along a regular closed curve C
is equal to the flux of the curl of that vector field through an arbitrary simply connected
(one-sided) surface S bounded by the curve C

 v dx =  curl v n dS . (2.81)
C S

2.4.10 Field types

There are several different types of fields depending on their properties, some of which are
listed below.
1) If s, v or T does not depend on time, such field is called stationary. Otherwise, the field is
non-stationary or transient.
2) If s, v or T depends only on two coordinates, the field is called two-dimensional, and if it
depends on one coordinate it is one-dimensional.
3) Vector field whose components are constant is called uniform.
4) Vector field v is potential (conservative) in a region V if it can be expressed as the gradient
of a scalar field, i.e. if at every point of V
v = grad s. (2.82)

Field s is called the potential of the field v. A potential field is also called non-circulatory as
curl v = curl (grad s) =   (s) = 0. (2.83)

This condition is a sufficient condition for the field to be potential.

5) A vector field v for which
div v = 0 . (2.84)

in every point of V, is called solenoidal (without source) in V.

19
 2. BASIC CONCEPTS AND DEFINITIONS

6) Vector field v which is potential and solenoidal at the same time is called Laplace or
harmonic field. For such a field:
v = grad s
2s
div v = div (grad s) =  (s) =  2 s = 0 or =0,
 x i x i

or in expanded form:
 2s  2s  2s
+ + = 0. (2.85)
x 2 y 2 z 2

The function s which satisfies this so called Laplace equation (2.85) is known as the Laplace
or harmonic function.

2.5 Elements of matrix algebra

It has been mentioned in the previous Section that vector and tensor components can be written
in the form of an array which is called a matrix. Also, systems of linear algebraic equations are
often associated with matrices. Matrices are very useful mathematical objects because they
enable an array of many numbers to be considered as a single object, denoting it by a single
symbol, and performing calculations with these symbols in a very compact form. Thus,
algebraic operations involving vectors and tensors, as well as the solution of a system of linear
algebraic equations can be conducted by employing the matrix algebra, a kind of 'mathematical
shorthand'.
In this Section some basics from matrix algebra, sufficient to follow this course, will be briefly
reviewed.

2.5.1 Definitions
1) A rectangular array of 'scalars' aij
a11 a12 ... a1n 
a 
 21 a22 ... a2n 
A = [aij]mxn = (2.86)
 ... ... ... ... 
 
am1 am2 ... amn 

is called an m x n matrix, or matrix of the format m x n, or simply a matrix. aij is the matrix
element situated in row i and column j, m is the number of rows, and n is the number of
columns. In general, number of rows m and columns n may not be the same.

2) n x 1 matrices are called column matrices or column vectors, and 1x n matrices are called
row matrices or row vectors. Their elements are usually written with one subscript only, e.g.

20
 2. BASIC CONCEPTS AND DEFINITIONS

 x1 
x 
 2
. 
x = [xi]nx1 =   , xT = [xi]1xn = [x1, x2, ..., xn]. (2.87)
. 
. 
 
 x n 

3) An n x n matrix is called a square matrix of order n. A square matrix A = [aij]nxn is:

– zero matrix if all its elements are zero, i.e. if aij = 0 for all i and j,
– diagonal matrix if all non-diagonal elements are zero, i.e. if aij = 0 for i  j ,
– identity or unit matrix if it is diagonal and if all diagonal elements are equal to 1, i.e.

I = [ij] = {1,0, if i= j
if i j (2.88)

where ij is known as the Kronecker delta.

– upper (lower) triangular if aij = 0 for i > j (i < j),
– strictly upper (lower) triangular if aij = 0 for i  j (i  j ).
4) The matrix AT = [aji]nxm, obtained from matrix A = [aij]mxn by interchanging rows and
columns is called the transpose of A (see (2.87)).
5) A (square) matrix is symmetric if it is equal to its transpose, i.e. AT = A, or aji = aij, and
skew-symmetric or anti-symmetric if AT = –A, or aji = –aij. Obviously, all diagonal elements
of a skew-symmetric matrix are zero.
6) Any matrix obtained by omitting some rows or/and columns from a given matrix A is called
a submatrix of A.
7) Although matrices do not have a numerical value, there is a number of 'scalars' associated
with matrices:
– For column and row matrices (2.87) one can define a number of norms, e.g. the
taxicab norm
n
||x||1 = | x i | , (2.89)
i =1
or the Euclidean norm
1
 n 2
||x||2 =   xi2  . (2.90)
 i =1 

– Similarly, for a square n x n matrix A, the taxicab norm is defined as

n n
||A||1 =   | aij | (2.91)
i =1 j =1

21
 2. BASIC CONCEPTS AND DEFINITIONS

and the Euclidean norm as

1
 n n 2
||A||2 =    aij2  . (2.92)
 
 i =1 j =1 

– Also, the trace of a square n x n matrix A is defined as the sum of its diagonal
elements
n
tr (A) =  aii , (2.93)
i =1

and the determinant of a square matrix A as

a11 a12 ... a1n
a21 a22 ... a2n
det (A) = . (2.94)
... ... ... ...
am1 am2 ... amn

8) The minor Dij of the element aij of a square matrix is the determinant of the submatrix
obtained by deleting the row and the column in which aij lies. The cofactor Aij of the
element aij is (–1)i+jDij.

9) The matrix A = [aij]mxn is said to have rank r if it contains at least one square submatrix of A
of the order r with a non-zero determinant, whereas the determinant of any submatrix of A
of the order r+1 or more is zero. The matrix A is non-singular or regular if and only if m = n
= r, i.e. if it is square and det (A)  0. Otherwise it is singular.

2.5.2 Matrix operations

Operations on matrices are defined in terms of operations on the matrix elements.

Equality
Two m x n matrices A = [aij]mxn and B = [bij]mxn are equal (A = B) if and only if aij = bij for all i
and j.
Addition and subtraction
The sum or difference of two m x n matrices A = [aij]mxn and B = [bij]mxn is the m x n matrix
A ± B = [aij ± bij]mxn. (2.95)

Scalar multiplication
The product of the m x n matrix A = [aij]mxn by the scalar s is the m x n matrix
sA = [saij]mxn. (2.96)

22
 2. BASIC CONCEPTS AND DEFINITIONS

Matrix Multiplication
The product of the m x n matrix A = [aij]mxn and the n x p matrix B = [bjk]nxp is the C = [cik]mxp
matrix
 n 
AB = [c ik ]m p =   aij b jk = [ai1b1k + ai 2b2k + ... + ainbnk ]m p . (2.97)
 j =1
  m p

It follows from the above definition that two matrices A and B can be multiplied in the order AB
if the number of columns in A is equal to the number of rows in B, i.e. if A and B are
conformable in the order AB. The existence of AB implies that of BA only if A and B are square
matrices. In general BA  AB.

Matrix Inversion
The inverse of the non-singular (square) n x n matrix A is the n x n matrix A–1 such that
A–1A = AA–1 = I, (2.98)
where I is the identity (unit) matrix. It can be shown that

[ ]n n = det1(A) [ A ji ]n n ,
1 1
A = aij (2.99)

where Aij is the cofactor of the element aij in the determinant det (A). The above definition
implies that an inverse of a singular matrix does not exist. If A–1 = AT, matrix A is said to be
orthogonal.
Useful relationships
It can be shown that the following relationships involving matrix operations hold
A + B = B + A, A + (B + C) = (A + B) + C
s (A + B) = sA + sB
A (B + C) = AB + AC
IB = B, CI = C (2.100)
A + 0 = A, 0A = 0, A + (–A) = A – A = 0
(AB)–1 = B–1A–1, (sA)–1 = s–1A–1, (A–1)–1 = A
(AB)T = BTAT, (AT)T = A
||A + B||  ||A|| + ||B||, ||sA|| = s ||A||, ||AB||  ||A|| ||B||,

where A, B and C are matrices of any format that allows the indicated operations, I and 0 are
identity and null matrices, respectively, and s is an arbitrary scalar.

Especially, for square matrices

IA = AI = A
tr (A + B) = tr (A) + tr (B), tr (sA) = s tr (A), tr (BA) = tr (AB)
det (AB) = det (BA) = det (A) det (B) (2.101)

23
 2. BASIC CONCEPTS AND DEFINITIONS

1 1
A= (A + A T )+ (A  A T ),
2 2
1 1
where, it should be noticed, (A + A T )is a symmetric and (A  A T )an asymmetric matrix.
2 2

2.5.3 System of linear algebraic equations

Definitions
A system of m linear algebraic equations with n unknowns x1, x2,..., xn

a11x1 + a12x2 + . . . + a1nxn = b1

a21x1 + a22x2 + . . . + a2nxn = b2
. . . . . . . . . . . . . (2.102a)
am1x1 + am2x2 + . . . + amnxn = bm

where aij are given numbers, which are called the coefficients of the system and bi are also
given numbers, can be written in a matrix form as
Ax = b, (2.102b)
where the coefficient matrix A = [aij]mxn, the solution vector x = [xi]nx1 and the right hand side
vector b = {bi}mx1 are
 x1  b1 
x  b 
a11 a12 ... a1n   2 2
a  .  . 
 21 a22 ... a2n 
A=
 ... , x =  , b =  . (2.103)
... ... ...  .  . 
  .  . 
am1 am2 ... amn     
 x n  bm 

The matrix
a11 a12 ... a1n b1 
a b2 
21 a22 ... a2n
B =  (2.104)
... ... ... ... ... 
 
am1 am2 ... amn bm 
is called the augmented matrix of the system (2.102).
Existence and general properties of solutions
System (2.102) may have no solutions, precisely one solution or infinitely many solutions.
Regarding the existence and the general properties of solutions of system (2.102), it can be
shown that:

1) System (2.102) has solutions if and only if the coefficient matrix A and the augmented
matrix B have the same rank r. If r = n, there is precisely one solution. If r < n, there are
infinitely many solutions which are obtained by determining r suitable unknowns (whose

24
 2. BASIC CONCEPTS AND DEFINITIONS

submatrix of coefficients has rank r) in terms of the remaining n – r unknowns, which can
take arbitrary values.
2) If system (2.102) is homogeneous (i.e. if b = 0), it always has the trivial solution x = 0. Non-
trivial solution exists only if rank of A is less than n.
Systems with square matrix
In a special case when the coefficient matrix A is a square n x n matrix, i.e. when the system
(2.102) has the same number of unknowns as equations, then

1) The system (2.102) has a unique solution

x = A–1b (2.105)

if and only if matrix A is non-singular, i.e. if det(A)  0 (which is the necessary and
sufficient condition for the existence of A–1),
2) The system of homogeneous equations has a non-trivial solutions if and only if the matrix A
is singular, i.e. if det (A) = 0.

Over-determined systems
Systems with more equations than unknowns (m > n) are said to be over-determined and, in
general, can be solved only approximately. One of the methods for solution of such systems is
the least square method. It consists in requiring that the sum of the squares of differences
between the right and the left hand side of each equation in (2.102a)
m
2
f(x1, x2,...,xn) = aij x j  bi =
2
 (aij x j  bi )2 (2.106)
i =1

has a minimum (note the use of the summation convention), i.e. that
m
f
= 2 aik (aij x j  bi ) = 0 (k = 1, 2, . . ., n), (2.107)
x k
i =1

which results in the following system of n linear algebraic equations with n unknowns
aikaijxj = aikbi (k = 1, 2, . . ., n) (2.108a)

(where i takes values from 1 to m, and j from 1 to n, m > n), or in the compact (easy to
remember) matrix form
ATA x = ATb. (2.108b)

Equation (2.108) is known as the normal equation of the system (2.102).

2.5.4 Eigenvectors and eigenvalues of a matrix

A column matrix (vector) x = [xi]nx1 is called the eigenvector of the square matrix A = [aij]nxn
if its elements are a non-trivial solution of the following system of homogeneous equations

25
 2. BASIC CONCEPTS AND DEFINITIONS

Ax = x, or (A – I)x = 0. (2.109)

A non-trivial solution will exist if and only if the determinant of the system vanishes
P() = det (A – I) = 0. (2.110)

The above non-linear algebraic equation is named the characteristic equation of A, and the
values of  which satisfy that equation are called the eignevalues of A. P() is called the
characteristic polynomial of A.

It can be shown that if matrix A is symmetric, all its eigenvalues are real, and if any two
eigenvalues are distinct, e.g. 1  2, then the two associated eigenvectors x1 and x2 are
orthogonal, i.e. x1T  x 2 = 0 .

Spectral radius and condition number

The set of eigenvalues of A is called the spectrum of A. The largest of the absolute values of the
eigenvalues of A, (A) = ||max, is called the spectral radius of A, while the ratio ||max/ ||min is
termed the condition number of A.

26
 3. STRESSES

Chapter 3

Stresses
3.1 Body and surface forces
Continuum mechanics is mainly concerned with forces acting on a continuum element and the
resulting motion. Those forces are conventionally classified into two categories: body (mass or
volume) and surface (traction) forces.

Body forces act on the total mass of the continuum element and are expressed as force per unit
mass of the element. These are, for example, gravitational force, inertial force, magnetic or
electro-magnetic force etc.
Surface forces are a consequence of the action of the immediate environment onto the
continuum element through their direct contact and are expressed as forces per unit area of the
element. These are, for example contact forces, forces of internal friction in a body etc. Surface
forces are also called stresses.

Different to particle mechanics, continuum mechanics is not concerned with forces themselves,
but with the density of their distribution in space (in a volume, or on a surface). Observe a
small volume of continuum no smaller than 3 (where  is of the order of magnitude of free
path of molecules) and call it a continuum element (Fig. 3.1).

3.1.1 Body forces

Thus, the concept of body force fb in point P of a continuum element understands the limit of
the ratio of the resultant force Fb acting on all points of that element and the mass of that
element m = V, when its volume V tends to zero (Fig. 3.1):
Fb Fb
fb = lim = lim . (3.1)
m0 m 
V 0 V
Consequently, the total body force acting on a given finite mass of continuum is
Fb =  fbdm =  fb dV . (3.2)
m V

27
 3. STRESSES

3.1.2 Surface forces

Analogously, the surface force in a point P is defined as
Fs
fs = lim , (3.3)
S 0 S
where Fs is the resultant of all forces acting on the surface S. According to the Cauchy
principal, this limit tends to a finite value dFs/dS while at the same time the moment of dFs
around the point P tends to zero. The total force acting on a finite surface is:
Fs =  fsdS . (3.4)
S

It is important to note an essential difference between the mass force fb and the surface force
fs. While the vector fb is a single-valued function of position vector and time, and consequently
forms a vector field, the vector fs has an infinite number of values in every point in space,
depending on the orientation of the surface on which it acts (Fig. 3.1), and, therefore, it does
not form any vector field.

a) b)

f sn
1
n1 f sn n2
2
S 1
P  S2 P

fb fb
V V

S S

Figure 3.1 Continuum element and body and surface forces

3.2 Stress tensor

In the previous sub-section it has been shown that, as opposed to body forces, the field theory
cannot be applied to surface forces, because vector fs does not form a field. However, it can be
shown that the surface force fs, acting on a surface with the unit outer normal vector n, and
henceforth denoted by fsn, can be expressed as the dot product of vector n and a (second order)
tensor  , which is a single-valued function of the position vector and time, i.e.
 =  (x,t ) , (3.5)

and therefore forms the field.

28
 3. STRESSES

Let’s consider a continuum element of a tetrahedral shape shown in Fig. 3.2. It is oriented in
such a way that three of its sides are in coordinate planes, while the position of the fourth side
is arbitrary and can be defined by its unit normal:
n = (nx, ny, nz) = (cosx, cosy, cosz), (3.6)

where:
nx = n  i = cos  x , n y = n  j = cos y , nz = n k = cos z . (3.7)

fs(-i)

n f sn
fs(-j)
y
S

fs(-k)
x
Figure 3.2 Cauchy tetrahedral

If fb is the resultant body force acting on the element and fsn, fs(-i), fs(-j), fs(-k) are surface forces
acting on the tetrahedral element sides, by applying the Newton second law of motion one
gets:
Va = Vfb + fsn S + fs (i) S x + fs ( j ) S y + fs (k) Sz , (3.8)

where a is the acceleration, V the volume and S, Sx, Sy, Sz are surfaces of the element
(Fig. 3.2). The term on the left hand side of eq. (3.8) and the first term on the right hand side
are proportional to the element volume. Being small quantities of the higher order they can be
neglected compared to other terms which are proportional to the surface, and eq.(3.8) becomes
fsn S = fsi Sx + fsj S y + fsk Sz . (3.9)

Since
S x Sy Sz
nx = cos  x = , n y = cos  = , nz = cos z = , (3.10)
S S S
then
fsn = n x fsi + n y fsj + n zfsk . (3.11)

By expressing the surface forces in terms of their components in Cartesian coordinate system:

29
 3. STRESSES

fsn = nx i + ny j + nz k ,

fsi = xx i + xy j + xz k ,
fsj = yx i + yy j + yz k , (3.12)
fsk = zx i + zy j + zz k ,

one gets:
nx = nx xx + ny yx + nz zx ,
ny = nx xy + ny yy + nz zy , (3.13)
nz = nx xz + ny yz + nz zz .
Expressions (3.13) show that components of fsn, which acts on an arbitrary surface, are
expressed as linear functions of the components of surface forces fsi, fsj, fsk, acting on
coordinate planes, i.e. in terms of nine quantities which can be regarded as components of a
tensor  with matrix:
 xx  xy  xz 
 
 =  yx  yy  yz  . (3.14)
 
 zx  zy  zz 

Tensor  is called (Cauchy) stress tensor. Consequently, eq. (3.11) can be written in the form:
fsn = n   , or fsn,i = njji . (3.15)

Thus, in every point of a continuum exists an infinite number of vectors fsn, whose direction
depends on the orientation n of the plane on which they act (e.g. Fig. 3.1), but only one stress
tensor  . The stress tensor, therefore, forms a (tensor) field. It is obvious that the components
of the stress tensor  , like the components of the body force vector fb, depend on the choice of
the coordinate system, but the tensor  , like the vector fb is an invariant (independent of the
coordinate system) physical quantity.
It will be shown later that stress tensor  is symmetric, i.e. that:
 =  T, or ij = ji . (3.16)

It is worth noticing that the first index in the expression for the surface force components (eqs.
(3.12) and (3.13)) denotes the direction of the normal to the surface on which the force acts,
while the second index denotes the direction of the axis to which the force is projected. For
example, zx denotes ‘x’ component of the force acting on the plane with normal in z direction.
The components with the same indices xx, yy, zz represent the projections of the force in the
direction of the surface normal and are called normal (or direct) stresses, while the projections
in the surface plane directions xy, xz, zy, … are called tangential (or shear) stresses (Fig.
3.3). Stresses can be positive and negative. The usual convention is that the stress is positive if
the outer normal and the force projection are both oriented in either positive or negative

30
 3. STRESSES

direction of the corresponding coordinate axes. Otherwise, the stress is negative. Positive
normal stresses cause tension while negative cause compression. All stresses in Fig. 3.3 are
positive.
z
zz k

zy j
zx i yz k
-yx i
- yy j xz k yy j
y
yx i
- yz k xy j
xxi

Figure 3.3 Normal and tangential stresses

In continuum mechanics, the stress components are often denoted with  for normal and  for
tangential stresses, i.e.
 xx  xy  xz 

 =  yx  yy  yz . (3.17)

 zx  zy  zz 

3.3 Principal stresses. Stress tensor’s invariants.

In general case surface force acting on the surface with normal n depends on the direction of n.
The angle between the surface force and the normal in a given point varies with the orientation
of the surface. It can be shown that a plane exists such that the surface force is in the direction
of the normal, i.e. a plane without tangential stresses. In fact, there are three planes orthogonal
to each other which satisfy this condition. Such plane is called a principal plane, and the
magnitude of the surface force (stress) on the principal plane is called principal stress.

If n is a unit normal to the principal plane (in the direction of the principal axis), then the
surface force acting on the principal plane is fp = n, where  is the principal stress. According
to (3.15) fp = n   . Since  is a symmetric tensor and n = In , then:
(   I) n = 0 or (ij - ij) nj = 0. (3.18)

31
 3. STRESSES

This system of three homogeneous linear algebraic equation has non-trivial roots (ni  0) if its
determinant is equal to zero, i.e. if:
det (ij - ij) = 0, (3.19)

which can be written in the following form:

3-I12-I2-I3 = 0 (3.20)

where:
 ) = ii ,
I1 = tr(
1
[ ]1
I2 = :  tr 2 () = [ ij  ij  ( ii ) 2 ] ,
2 2
(3.21)
1
I3 = det( ) = ijk pqr ip jq kr ,
6
are the first, second and third invariants of the stress tensor.
If 1, 2 and 3 are the roots of eq. (3.20), and therefore principal stresses, then eq. (3.20) can
be written in the form:
( -1) ( -2) ( -3)=0, (3.22)

and it follows that:

I1 = 1 + 2 +3 ,
I2 = -(12 + 23 + 31) , (3.23)
I3 = 123.

For each i the system (3.18) has a non-trivial solution which, using the condition nini = 1,
gives the components of one principal direction.
If principal directions are chosen as coordinate axes, stress tensor becomes diagonal (cannonic
form):
 1 0 0 
 
 = 0  2 0  . (3.24)
0 0  3 

3.4 Spherical and deviatoric stresses

Some materials are more sensitive to stress components which cause the change of shape than
to those causing the change of volume. Examples include nearly all fluids and many solids
experiencing plastic deformations. In such cases it is desirable to define the stress tensor that
would be the same for all stress states which differ only by the pressure to which a material is
exposed.

This is achieved by defining a mean normal stress as:

32
 3. STRESSES

1 I
p = -p= ii = 1 ., (3.25)
3 3
where conventionally pressure p has a negative sign since the stresses of compression
(pressure) are negative. Then, the stress tensor can be divided into spherical and deviatoric
parts, i.e. into the sum of hydrostatic or spherical stress which causes the change of volume,
and deviatoric stress which causes the change of shape (in isotropic materials):
 = pI +  d = [pij + (ij – pij)] iiij . (3.26)

The principal directions of  and  d are the same since they represent the normals to the
planes with zero tangential stresses, and their principal values are:
i =  p +  di (i = 1, 2, 3). (3.27)

33
 4. DEFORMATION AND FLOW

Chapter 4

Deformation and flow

4.1 Introduction
In continuum mechanics, materials are generally classified as solids and fluids, depending on
their behaviour when subjected to loading (mechanical and/or thermal). A solid body deforms
until its internal forces are in equilibrium with externally applied loads, while a fluid body flows
under the action of shear stresses. However, a fluid body also deforms during the flow, and a
solid body is said to flow if the deformations are large (hyper-elasticity, plasticity, …). There
are also materials with characteristics of both solids and fluids (such as visco-elastic or visco-
plastic materials) and cannot be easily classified.

4.1.1 Point and particle. Body and configuration

In the kinematics of continua, one must distinguish between the point in the space and the
material point. A particular position in space is called a point while a small part of a continuum
is called a particle. The (material) body is then defined as the system of particles, and the
domain in space occupied by the body at a given time is called the current configuration.

Deformation is the change in the dimensions and shape of a body between the initial (reference,
undeformed) and final (deformed) configurations. Generally, only initial and final
configurations are considered, without paying attention to intermediate stages of deformations.

Flow is defined as a continuous motion of a continuum and is represented by a velocity field.

4.1.2 Configurations
Let’s assume that at a reference time t = t 0 , which is usually taken as zero and is called the
initial time, a material body B occupies the part (region) of the space R0. Particles of body B
continuously fill the space R0, and there is a one-to-one correspondence with particles of body B
and points in R0. The position of each particle P0 of the body can then be defined with the
position vector  (Fig. 4.1), which in a fixed coordinate system O0123 is given as:
 = 1i1(0) + 2i (0) (0)
2 + 3i 3 . (4.1)

34
 4. DEFORMATION AND FLOW

In the next time instant t the body B continuously occupies the space R, and all particles of B
uniquely correspond to points in R. Then, a new position of the particle P0, denoted as P in the
new configuration R, can be specified in a coordinate system Ox1x2x3 by a position vector:
x = x1i1 + x 2i 2 + x 3i 3 . (4.2)

3
x3

P0
P

(0)
x
i3 i3
i2
O0 O x2
i1
i1
(0)
i2
(0)
2

1
x1
Initial (reference, undeformed) Current (deformed)
configuration configuration

Figure 4.1 Initial and current configurations

4.1.3 Material and spatial description

Since the particle P0 in the reference configuration R0 with position vector  takes the position P
in the current configuration R with the vector position x, the motion of body B can be described
by:
x = x(,t) or xi = xi ( k ,t) . (4.3)

This expression gives the current location x of the particle which was at point  at time t = 0.
Equation (4.3) can be interpreted as a mapping of the initial configuration into the current
configuration. If the Jacobian of this transformation
 x 
J = det i
 0 , (4.4)
  j

then the inverse transformation exists such that:

 = (x,t) or i = i (x k ,t) . (4.5)

A typical particle which moves according to (4.3) must have its reference (initial) position  in
R0 and, therefore, point  is its reference point. Regardless of the particle motion, it can only
change its current position and it will always have the same reference position. Coordinates i
can therefore be regarded as fixed and are called material or Lagrangian coordinates. The

35
 4. DEFORMATION AND FLOW

system of all particles at t = t 0 , i.e. region R0, is called the reference (initial or undeformed)
configuration.

On the other hand, the position vector x denotes the point in space which can be ‘filled’ with
different particles at different instants of time. Coordinates xi are called spatial or Eulerian
coordinates, and the region R occupied by the body at time t is called the current or deformed
configuration.
Continuum motion can be described by using either material i or spatial xi coordinates. The
current position x of a particle at any time t can be found if its reference position  is known, or
the reference position  can be found if the current position x is known. If the motion is
desribed using eq. (4.3) such a description is called a material or Lagrangian description, and
the description by (4.5) is called spatial or Eulearian.

When using the Lagrangian description, a particle is followed in time. Eulerian description
gives information of the motion in a given spatial point of different particles which ‘fill’ that
point at different times. The choice of the description will depend on the problem considered,
but usually the Lagrangian description is used for solid body motion (and deformation), and the
Eulerian for fluid body motion.

4.1.4 Displacement vector

The vector connecting the points P0 and P (Fig. 4.1) is called the displacement vector and can
be expressed in either material or spatial coordinates:
r (0) (0)
u = P0 P = ui i i = ui ii , (4.6)
r
or in terms of x,  , and vector O0O :
r
u = x   + O0O . (4.7)

In continuum mechanics, the coordinate systems O0123 and Ox1x2x3 are normally chosen to
coincide to each other, and then:
u= x . (4.8)

In the Lagrangian description, displacement u is represented as a function of  and t:

u(, t) = x(, t)   , (4.9)

and in the Eulerian description as function of x and t:

u(x, t) = x  (x, t) . (4.10)

These expressions have different physical interpretation. Expression (4.9) represents the
displacement at time t of the particle with the position vector  , while expression (4.10)
represents the displacement of any particle which is in the position x at time t.

36
 4. DEFORMATION AND FLOW

4.2 Deformation

4.2.1 Measure of deformation

When considering a solid body under the action of forces, often the main aim is to find its
deformation. This requires some kind of deformation measure to be defined. If, for example,
one considers a 1D rod of the initial length L0 stretched to a length L, then the total change of its
length L = L - L0 can be used as the measure of the deformation. It is obvious that under the
same load higher L will correspond to a longer rod. Therefore, a relative change in the
dimension is used as the measure of deformation, such as:
L L  L0 L
= , n = = ,
L0 L0 L0
L
dL  L
t =  L
= ln = ln( ) = ln(1+  n ),
 L0 
(4.11)
L0
L L0 L
 = =
L L
where  is the lengths ratio (stretch), and n (engineering strain), t (true strain) and ’ are
measures of relative extensions. Also, the following measures are sometimes more convenient:
L2  L20 L2  L20
e= , e = . (4.12)
2L2 2L20

For small deformations all of these measures are practically equivalent. The choice of a
particular measure of deformation will depend on the relationship between the stress and the
deformation, i.e. on the constitutive relationship.

4.2.2 Gradients

Deformation gradients

Let’s observe an infinitesimal distance between two particles of a continuum body subjected to
deformation (Fig. 4.2). Then, for example:
x
dx = F  d or dx i = i d j , (4.13)
 j

and:
i
d = H dx or di = dx , (4.14)
x j j

where

37
 4. DEFORMATION AND FLOW

xi
F = (grad x )T or Fij = , (4.15)
 j

is the material deformation gradient, while

T 
H = (grad x ) or Hij = i , (4.16)
x j

is the spatial deformation gradient. It is obvious that:

F H = H F = I , (4.17)
and therefore , H = F-1 is inverse of the tensor F. It could also be noticed that:
1
det(F) = =J . (4.18)
det(H)

Q
x3 3 u+du
Qo

+d  dx
d
u
 Po P
x
x2

x1 2
1

Figure 4.2 Deformed and undeformed configurations

Displacement gradients
If the relative displacement du = dx - d is taken as the measure of the deformation, one gets:
du = dx  d = F  d - I d = (F  I) d = J d , (4.19)

or
du = dx  d = I dx - H dx = (I  H) dx = K  dx , (4.20)

where:
u x
J = (grad u )T = [grad (x  ) ]T = F  I or Jij = i = i   ij , (4.21)
 j  j

is the material displacement gradient, and

T T u 
K = (grad x u ) = [grad x (x   ) ] = I  H or Kij = i = ij  i , (4.22)
x j x j

38
 4. DEFORMATION AND FLOW

is the spatial displacement gradient.

4.2.3 Deformation tensors

If either F or J (or indeed H or K) are known for a body at every time instant, then the local
deformation is described for all particles of the body. However, F and J (as well as H and K)
include not only deformation but also rigid body rotation which does not contribute to the
change of shape nor to the stresses in the body (as will be shown later). Therefore it is
convenient to define the deformation measures which ignore rotations. Thus, instead of working
T T
with vector quantities dx and d one should use their invariants dx  dx or d  d.

Green and Cauchy deformation tensors

According to (4.13)
2 T T T T
(dx) = dx dx = (F  d)  (F d) = (d F ) (F  d) , (4.23)

or
x k x k
(dx) 2 = dT G  d or (dx)2 = d d , (4.24)
i  j i j

where
T x k x k
G = F F or Gij = (4.25)
i  j

is the Green deformation tensor. Similarly according to (4.14)

2 T T T T
(d) = d d = (H dx)  (H dx) = (dx H ) (H  dx) , (4.26)

or
k k
(d) 2 = dxT C  dx or (d)2 = dx dx , (4.27)
xi x j i j

where
T k k
C = H H or Cij = (4.28)
xi x j

is the Cauchy deformation tensor. Similar to tensors F and H, tensors G and C are inverse of
each other, i.e. C = G-1.

Lagrange and Euler tensors of finite strain

The difference between (dx)2 and (d )2 can be taken as the local measure of deformation
between initial and final configurations:
2 2 T T T T
(dx)  (d) = d G  d  d  I d = d  (G  I)  d = d 2L d , (4.29)

39
 4. DEFORMATION AND FLOW

where:

1 1 1  x x

L= (G  I) = (FT  F  I) or Lij = k k  ij  , (4.30)

2 2 2  i  j 

is Lagrange or Green finite strain tensor. Similarly

2 2 T T T T
(dx)  (d) = dx I  dx  dx  C dx = dx  (I  C)  dx = dx 2E dx , (4.31)

where:

1 1 1  

E= (I  C) = (I  HT  H) or Eij =  ij  k k  , (4.32)
2 2 2 x i x j 

is Euler or Almansi finite strain tensor.

According to (4.21), F = J + I, and

L=
1
2 ( ) (
1
( J + I)T  (J + I)  I = JT + J + J  JT , or
2 )
L=
1

2
( T
) ( T
)
grad u + grad u + gradu  grad u , or

(4.33)

1  u i u j u k u k 
Lij = + +
2   j i i  j 
.

Similarly, according to (4.22), H = I – K, and

E=
1
2 ( 1
) (
I  (I  K )T  (I  K ) = KT + K  K K T , or
2 )
E=
1
2
( [ T T
]
grad x u) + grad x u  grad x u  (grad x u) , or (4.34)

1  ui u j uk uk 
E ij = + 
2  x j xi x i x j 
.

Factor 1/2 is introduced so that L and E would reduce to the corresponding small strain tensors
after higher order terms are neglected. It can be noticed that in cases where there is no
deformations (du = 0, dx = d ), L and E become null tensors while G and C are unity tensors.

4.2.4 Small deformations

If the deformations are very small (infinitesimal), the above relations can be simplified. In cases
where det(J)<<1 or det(K)<<1, which is equivalent to du << dx , or du << d , or dx  d,
T T T
then the terms J J or K K , being the small higher order terms compared to J + J or
K T + K , can be neglected. Also, the difference between the undeformed and deformed

40
 4. DEFORMATION AND FLOW

configurations can be ignored (dx  d), and L and E reduce to the well known Cauchy small
(or infinitesimal) strain tensor:

Ls = E s =
1
2 [
(gradu)T + gradu , or ]
1  u u j  1  ui u j 
lij = ij =  i + = +
2   j i
2  x j xi

, (4.35)

where the grad operator can be applied in either Lagrangian or Eulerian coordinates, and lower
case letters l and  are used to denote small strains.

Expressions (4.35) can also be obtained in the following way. If u = u(x,t) is expanded in the
vicinity of P0 by using the Taylor series expansion and neglecting the higher order terms, one
gets:
 u 
u = u0 + du = u0 + (grad u)0  dx or ui = ui 0 +  i  dx j .
T
(4.36)
 x j  0

The tensor (grad u)0 can be expressed as the sum of its symmetric and anti-symmetric parts
T

giving:

u = u0 +
1
2
[ 0
1
2
] [
(grad u)T + grad u dx + (grad u)T  grad u dx
0 ,
] (4.37)
= u0 + E s0  dx + 0 dx

where E s0 is Cauchy small strain tensor and

=
1
2[(grad u)T grad u ] 1  u
or ij =  i
u j 
2  x j xi

, (4.38)

is linear rotation tensor, and:

K = (grad u) = E +  or Kij = ij + ij .
sT
(4.39)

If  is expressed in term of its rotation vector

1 1 u 1
 = curl u or  i = eijk k = eijk kj , (4.40)
2 2 x j 2

then expression (4.37) can be written as:

s
u = u0 +   dx + E0 dx . (4.41)
Rigid body Deformation
motion

It can be seen that the motion of a body can be decomposed into two parts: rigid body motion
consisting of translation and rotation, and deformation. The first part obviously does not cause
stresses and is therefore neglected in continuum mechanics.

41
 4. DEFORMATION AND FLOW

Diagonal elements of the small strain tensor represent the change in length per unit original
length in the direction of the coordinate axes, and are called normal strains. Off-diagonal terms
represent one-half the angle between two line elements originally at right angle to each other,
and are called shear strains. The concept of small (infinitesimal) deformations is used in the
analysis of elastic bodies. It is also useful as a first step in analysing of large deformations.

4.2.5 Principal strains. Invariants of strain tensors.

Lagrangian and Eulerian small strain tensors are symmetric tensors of the second order, and
their principal values, principal directions and principal planes can be obtained in the way
described in Section 2.4.7 and 3.3. For example, characteristic equation for Lagrangian small
strain tensor
det(Ls-l I)=0 or det (lij - lij) = 0, (4.42)

can be written in the following form:

l3-I1 l 2-I2 l -I3 = 0 (4.43)
where:
I1 = tr(Ls) = lii = l1 + l 2 + l 3,

2[
1 s s 2 s
] 1
2
2
I2 = L :L  tr (L ) = [lij lij  (lii ) ] = ( l1l2 + l2l3 + l3l1) , (4.44)
1
I3 = det(Ls) = eijke pqr lip l jqlkr = l1l2l3 ,
6
are the invariants , and the roots of (4.43) l1, l2, l3 are principal values of the Lagrangian small
strain tensor. Principal directions and planes correspond to these principal values. Analogue
relationships can be obtained for Eulerian small strain tensor Es.
It is worth noticing that small deformation tensors L, E, G and C are symmetric and therefore
have three real principal values with corresponding three (orthogonal) principal directions.
Since L and G, and E and C differ by the constant, their corresponding principal directions
coincide.

4.2.6 Polar decomposition of the finite deformation tensor

Decomposition of the finite deformation tensor F ( Fij = x i /  j ) into the sum of its symmetric
and antisymmetric parts represents the decomposition of F into the sum of pure rotation and
pure stretch. In order to eliminate rigid body rotation which does not affect the stress state in the
body, it is sometime convenient to use polar decomposition of F into the product of two tensors,
one representing the rigid body rotation and the other representing a stretch and is symmetric. If
R is an orthogonal rotation tensor, which rotates principal axes of the Green deformation tensor
G in point  onto the principal axes of Cauchy deformation tensor C in x, then two symmetric,
positive definite tensors U and T exist such that:

42
 4. DEFORMATION AND FLOW

xi
F = R  U = T R or Fij = = RikU kj = Tik Rkj , (4.45)
 j
where U is the right and T is the left stretch tensors. This decomposition can be interpreted
using expression (4.13) and
dx = F  d = R U  d , (4.46)

as the body stretch by U followed with the rigid body rotation by R (Fig. 4.3), or as the rigid
rotation by R followed by the stretch T.
F

d dx

U R
,
dx

Figure 4.3 Interpretation of the polar decomposition

Polar decomposition is theoretically very important, but the stretch tensors are not practical
since their components are complex functions of deformation gradients. From the practical
calculations, components of the deformation tensors G and C are more useful.

4.2.7 Dilatational and deviatoric components of deformation tensor

If the mean volumetric dilatation is defined as:
I 1
v = 1 = lii , (4.47)
3 3
then Lagrangian small deformation tensor can be decomposed into the sum of its spherical
(dilatational) and deviatoric tensors in the following way:

[ ( )]
Ls = Ls,s + Ls,d =  v I + Ls,d =  v ij + lij   v ij i ii j . (4.48)
The spherical component represents the change in volume and it can be shown that:
V0
 l1 + l2 + l3 = I1 = 3 v , (4.49)
V0
where V0 is the initial volume. The deviatoric component represents the change in the shape of a
body.

43
 4. DEFORMATION AND FLOW

4.2.8 Compatibility conditions

Compatibility conditions require that kinematic variables which describe the motion
(displacements, displacements’ gradients,…) are compatible with the assumption of single-
valued, continuous motion x = x( , t). Physically, this means that a continuum particle cannot
be in two different places at a given time, nor that two particles can be in the same position at
the same time (no holes nor overlapping in the continuum body).

Compatibility is automatically satisfied if a problem is mathematically defined and solved in

terms of displacements functions, which are continuous. However, it is sometimes more
convenient to formulate the problem in terms of strains (such as in problems of linear elasticity).
In this case, continuity of displacements does not mean the continuity of strains and additional
compatibility conditions are required.

If, for example, the components of Cauchy small strain tensor are expressed explicitly in terms
of coordinates, then six independent equations (4.35)
1  u u j 
ij =  i +
2  x j xi 
(4.50)

can be regarded as six partial differential equations with three unknown displacement function
ui. This system is obviously over-determined and, in general, its solution does not exist for
arbitrary ij, or it is not compatible, and additional conditions are needed.

By differentiating expressions (4.50), one gets:

1  3ui  uj 
2 3
 ij
=  + , (4.51)
x k x l 2  x jx k x l xi xk x l 

and by replacing indices:

2 kl 1  3uk 3ul 

=  + ,
x ix j 2  x l xi x j x k xix j 

2 jl 1  uj 3ul 
3
=  + , (4.52)
x ix k 2  x l xi x k x j xix k 

2ik 1  3ui  3uk 

=  + .
x j x l 2  x k x jx l xi x jx l 

It follows from the symmetry of Es and the continuity of ij that:

2 2
 ij 2 kl  2 ik   jl
+   = 0. (4.53)
x k x l x ix j x j x l xi x k

44
 4. DEFORMATION AND FLOW

This is the necessary and satisfactory compatibility conditions, known as St. Venant
compatibility equations. There are 81 equations but only 6 are independent (due to symmetry of
Es):
211 2 22 212
+ = 2 ,
x 22 x12 x1x 2
2 22  2 33  2 23
+ =2 ,
x32 x 22 x2x 3
2 2
  33  11  2 31
+ = 2 ,
x12 x 32 x3x1
(4.54)
   23  31 12  211
 + +
= ,
x1  x1 x 2 x3 x2x 3
2
   23  31 12    22
  +
= ,
x 2  x1 x 2 x3 x3x1
   23  31 12  233
 + +
= .
x 3  x1 x 2 x3 x1x 2

In case of finite strains, the compatibility conditions are considerably more complex.

4.3 Motion and flow

Motion and flow describe instantaneous or continuous changes of the configuration. In solids,
flow is used to define the permanent deformation of solids, while in fluids it means continuous
motion.

The Lagrangian description follows the continuum particle and is conventionally used in solid
mechanics where the deformation is of particular interest. On the other hand, the Eulerian
approach is used in fluid mechanics where the velocity field is the main information required.
While the Lagrangian approach focuses on continuum particles (defined by their positions at t =
0), the Eulerian approach is concerned with a particular region of space (control volume).

Since the functions (4.3) and (4.5) are inverse of each other, any physical (scalar, vector or
tensor) continuum quantity F expressed for a continuum particle (in Lagrangian or material
description) with a function
F = F1(,t) = F1(i ,t) , (4.55)

can also be expressed in reference to a particular location in space occupied by the particle
(Eulerian or spatial description) with a function
F = F1((x,t),t) = F2 (x,t) = F2 (xi ,t), (4.56)

45
 4. DEFORMATION AND FLOW

where indices 1 and 2 in F1 and F2 are employed to demonstrate that a function in Eulerian
coordinates is not in general the same as in Lagrangian coordinates.

4.3.1 Rate of deformation and material derivative

The rate (speed) at which a certain continuum property changes in time, i.e. the rate of change
that would be observed by the observer travelling with the particle under consideration, is called
material derivative. Time derivatives with material coordinate  held constant are called
material time derivatives. It is sometimes denoted by d/dt, but more often by D/Dt in order to
underline that the derivative is material derivative.

The instantaneous position vector of a given particle is a function of time, and its material
derivative is the velocity of that particle and, at the same time, represents the velocity of
continuum in a given point of space, i.e.
dx Dx
= = v(x,t). (4.57)
dt Dt
This is an important relation, especially in fluid mechanics, where the function  = (x,t) , i.e.
the initial positions of fluid elements passing through a given domain of space (control volume)
is neither of interest nor important.
In Lagrangian coordinates
DF DF1 (,t) F1 (,t)
= = , (4.58)
Dt Dt t
while in Eulerian coordinates (notice that x is a function of t with material derivative Dx/Dt)
DF DF2 (x,t) F2 Dx
= = + gradF2  , (4.59)
Dt Dt t Dt
or by using (4.57)

DF F
= + v(x,t) gradF . (4.60)
Dt t
Local Material
change change

The first term on the right hand side of expression (4.60) represents the temporal rate of change
in a given point of space (local change) and is also called spatial time derivative. Time
derivatives with spatial coordinate x held fixed are called spatial time derivatives. The second
term is a consequence of the fact that the continuum changes the position in space (material
change) and is called convective or transport term.

46
 4. DEFORMATION AND FLOW

Based on (4.60), one can define operators for material derivative in Lagrangian coordinates:
D(...) (...)
= . (4.61)
Dt t
or in Eulerian coordinates:
D(...) (...)
= + v grad(...) . (4.62)
Dt t

4.3.2 Velocity and acceleration

One definition of velocity is given by (4.57). Alternatively, if x is given as  + u (eq. (4.8)), the
velocity can be expressed via displacement vector:
Dx D( + u) Du
v= = = . (4.63)
Dt Dt Dt
In Lagrangian coordinates u = u( ,t) so
x(,t) u(, t)
v(, t) = = , (4.64)
t t
while in Eulerian coordinates u = u(x,t) and the velocity is given with the following implicit
formula:
u(x, t)
v(x, t) = + v(x, t) grad u(x, t) . (4.65)
t
This expression is usually not employed since the so called instantaneous velocity
v = v(x, t), (4.66)

is used as dependent variable (in fluid mechanics) and not x or u.

Material derivative of velocity is acceleration. According to (4.61) and (4.64), the acceleration
in Lagrangian coordinates is given as:
Dv v  2x  2u
a(, t) = = = = , (4.67)
Dt t t 2 t 2

while in Eulerian coordinates according to (4.62):

Dv v
a(x, t) = = + v grad v. (4.68)
Dt t

In relation to Lagrangian description, one can define particle trajectory as the system of all
spatial points passed through by the particle. With reference to Eulerian description, a
streamline is defined as the line of the velocity vector field, i.e. at every time instant velocity
vectors are in the direction of tangent to the line. In a given time instant, a velocity vector is
tangent to both trajectory and streamline. The trajectory equation is obtained by solving the
following vector equation:

47
 4. DEFORMATION AND FLOW

dx
= v(x,t). (4.69)
dt

with the initial condition x =  at t = t0, while the vector equation for the streamline at t = t0 is:
dx = v(x,t 0 ). (4.70)

where  is a parameter. If the velocity is not function of time, the corresponding streamline and
trajectory coincide.

4.3.3 Material derivative of a volume integral. Reynolds transport

theorem.
The material time derivative of an integral is the rate of change of that integral on a material
domain. If F is an arbitrary single-valued function of spatial coordinates and time with
continuous derivatives, then:

 FdV =  F (x,t)dV , (4.71)

V V

where V = V(t) is the material volume (i.e. the volume which moves with the continuum and no
mass flux occurs across the volume boundary), is also well defined function of space and time.
By substituting variables in (4.71) and using (4.3), one gets:

 F(x,t)dV =  F[x(,t),t]JdV0 , (4.72)

V V0

where J = dV/dV0 is the Jacobian of the transformation (4.3) defined by expression (4.4), and V0
= V(t0) is the material volume at t = t0. Based on (4.72), the material derivative of the integral
(4.71) can be expressed as:
D D
Dt  F (x,t)dV =
Dt  F[x(,t),t]JdV0 . (4.73)
V V0

Noticing that V0 does not depend on time, which means that the order of differentiation and
integration can be swapped, and using Euler expansion formula:
DJ
= Jdivv . (4.74)
Dt
one gets:
D  DF DJ   DF 
Dt    Dt   Dt
FdV = J +F dV = + Fdivv JdV0 . (4.75)
Dt  0 
V V0 V0

48
 4. DEFORMATION AND FLOW

Returning to Eulerian coordinates and using (4.60), expression (4.75) becomes:

D  F  F 
Dt  FdV =  
t
+ v gradF + Fdivv dV =   + div(Fv) dV .

t
(4.76)
V V V

Finally, by employing Gauss theorem (2.80) one gets:

D F
Dt  FdV = 
t
dV +  Fv ndS , (4.77)
V V S

where S is the area bounding the material volume V. The equation (4.77) represents the
Reynolds transport theorem which states that the rate of change of integral of function F over
the material volume, is equal to the sum of the rate of change of that integral over the fixed
volume in space which coincide with the material volume at a given time, and the flux of F
through the bounding surface S.

4.3.4 Rate of deformation

To obtain the relation between the rate of deformation and the velocity of a continuum particle,
consider two close points P0 and P, with position vectors x0 and x = x0 + dx, and velocities v0
and v, respectively (Fig. 4.4).

x3 P
v

dx dv

x
Po
vo
xo

x1 x2

Figure 4.4 Velocities of the continuum particles

If the function v(x, t) is expanded in the vicinity of P0 by using the Taylor series expansion and
higher order terms are neglected, one gets:
 v 
v = v 0 + dv = v0 + (grad v)0  dx or vi = vi 0 +  i dx j ,
T
(4.78)
 x j  0
where
v i
Y = (grad v)
T
or Yij = , (4.79)
x j
is the spatial velocity gradient.

49
 4. DEFORMATION AND FLOW

As before, the tensor (grad v )T0 can be decomposed into the sum of its symmetric and anti-
symmetric parts giving:

v = v0 +
1
2
[ 0
1
2
] [
(grad v)T + grad v dx + (grad v)T  grad v dx
0 ,
] (4.80)
= v 0 + D0 dx + V0  dx

where

D=
1
2[(grad v)T + grad v ] 1  v
or Dij =
i +
v j 
2  x j xi 
, (4.81)

is the deformation rate tensor, and

V=
1
2[( grad v)  grad v
T
] or Vij =
1  vi v j 

2  x j xi 
, (4.82)

is the vorticity or spin tensor, and:

Y = (grad x v) = D + V or Yij = Dij + Vij .

T
(4.83)

If V is expressed in term of its vorticity vector

1 1 v 1
˙ = curl v or ˙ i = ijk k = ijkVkj , (4.84)
2 2 x j 2

then expression (4.78) can be written as:

v = v 0 + ˙  dx + D0  dx. (4.85)
Rigid body Deformation
motion

The rate of deformation tensor D is not limited to small deformations although it has the same
form as the small strain tensor Es (or Ls). It can be seen that, analogously to deformation, the
velocity of a body can be decomposed into two parts: the first corresponding to the velocity of
rigid body motion (translation and rotation), and the second corresponding to the rate of
deformation. The first part does not cause stresses and is normally neglected in continuum
mechanics.

It is easily shown that:

[ ]
s
DL DE s 1
= (grad v) + grad v = D ,
T
= (4.86)
Dt Dt 2

and

D 1
[
= (grad v)  grad v = V .
Dt 2
T
] (4.87)

Expression (4.86) can be written as:

50
 4. DEFORMATION AND FLOW

DEs = DDt or D ij = Dij Dt . (4.88)

where components Dij are known as the true (natural) strain increments.

4.3.5 Principal values and invariants of the rate of deformation

tensor.
The rate of deformation tensor D is symmetric tensor of the second order, and its principal
values, principal directions and principal planes can be obtained in the way described in Section
4.2.5. The characteristic equation
det(D-DI) = 0 or det (Dij - Dij) = 0, (4.89)

can be written in the following form:

D3-I1 D2-I2 D -I3 = 0, (4.90)
where:
I1 = tr(D) = Dii = D1 + D 2 + D3,
[
1
2
2
] 1
2
2
I2 = D:D  tr (D) = [Dij Dij  (Dii ) ] = (D1D2 + D2 D3 + D3D1 ) , (4.91)
1
I3 = det(D) = ijk pqr Dip DjqDkr = D1D2D3 ,
6
are the invariants , and the roots of (4.90) D1, D2, D3 are principal values of the rate of
deformation tensor. Principal directions and planes correspond to these principal values.

4.3.6 Dilatational and deviatoric components of the rate of

deformation tensor
An important parameter in the analysis of continuum motion is the rate of volumetric dilatation,
which is defined as the rate of volume change per unit volume. By employing (4.74) one gets:
1 DV 1 D( JV0 ) 1 DJ
= = = divv = tr(D) = I1 , (4.92)
V Dt JV0 Dt J Dt

showing that the rate of volumetric dilatation is equal to the sum of diagonal elements of the
rate of deformation tensor.
If the mean rate of volumetric dilatation is defined as:
I 1
 = 1 = Dii , (4.93)
3 3
then the deformation rate tensor can be decomposed into the sum of the spherical (dilatational)
part, which corresponds to the rate of the volume change, and the deviatoric part which
represents the rate of change of shape:
s d d
(
D = D + D = I + D or Dij = ij + Dij  ij . ) (4.94)

51
 5. FUNDAMENTAL LAWS

Chapter 5

Fundamental laws of continuum mechanics

5.1 Introduction
In the previous two chapters the forces (body and surface) and the motion (deformation and
flow) were introduced independently of each other. In this chapter, the forces will be linked
to the resulting motion. For that purpose, the following fundamental laws of physics will be
used:
1) Mass conservation,
2) Conservation of linear momentum (first Euler’s law, second Newton’s law),
3) Conservation of angular momentum or moment of momentum law (second Euler’s law),
4) Energy conservation (first law of thermodynamics),
5) Entropy production (second law of thermodynamics).

Apart from the fundamental laws, which are independent of the type of continuum under
consideration, additional laws or constitutive relationships will also be used. These are
specific for each continuum, such as Hooke’s elasticity law, Stoke’s viscosity law, Fourier’s
law of heat transfer, etc.

5.1.1 General transport equation

A material domain or a system is defined in the previous chapter as an assemblage of
material particles. A material domain moves with the material, so that material points on the
boundary remain on the boundary and no mass flux occurs across the boundary. A control
volume is defined as a fixed region of space occupied by the material at a given time. It is
easy to notice that the system is related to the Lagrangian description of motion while the
control volume is related to the Eulerian description. The two are linked via Reynold’s
transport theorem (4.77). In this chapter, equations will be given for both the system and the
control volume in the original integral form and also in differential form.

52
 5. FUNDAMENTAL LAWS

For a system:
If  denotes a global property of a system such as mass, or momentum and  =  (x, t)
represents the density of that property, then the relationship between them is given by a
material integral (over a given material domain instead of a given spatial domain) as:
H=  dV , (5.1)
V

where  is the mass density and V is the system volume (material domain). All fundamental
laws are concerned with the rate of change of  due to the interaction between the system
and its environment. These changes can be proportional to the mass of the system or to the
surface bounding the system as:
DH
Dt
=   m dV +   S
 ndS . (5.2)
V S

This expression can be regarded as the general conservation law or the general transport
equation. By using Gauss theorem (2.80) one gets:

 (  m + div  ) dV ,
DH
= (5.3)
Dt S
V

and due to (5.1) [and (5.20)], the differential form of the general conservation law becomes:
D(  )
=  m + div  S . (5.4)
Dt
Volumetric (body, mass) term m and surface term  s will be given in the next section for each
fundamental law.

For a control volume:

Solid mechanics normally uses the system since it is easily defined for a solid body. On the
other hand, it is not easy to define the system for a fluid flow, and it is often more convenient
to consider a region of space through which the fluid flows, so called control volume (CV).
The size and shape of the control volume are arbitrary, and are usually chosen so that control
surface (CS) partly coincides with flow boundaries. The control volume is fixed to a chosen
coordinate system, which is not necessarily fixed.
Taking that F =  in the Reynold’s transport theorem (4.77), equation (5.2) at the moment
when the system (material domain) coincides with the control volume becomes:
 ( )
 dV +  ( v n)dS =   m dV +   S ndS , (5.5)
t
V S V S

where V and S represent the volume and the bounding surface of the control volume,
respectively. This is the general conservation law for the control volume in the integral form.

53
 5. FUNDAMENTAL LAWS

By applying the Gauss theorem again, one gets:

  (  )

  t 
( )
+ div (v) dV =   m + div  S dV . (5.6)
V V

Since V is arbitrary, then:

(  )
+ div (v) =  m + div  S . (5.7)
t
which is a differential form of the general conservation law for the control volume.

5.2 Fundamental laws of continuum mechanics in various forms

Expressions for fundamentals laws of continuum mechanics can now be derived from the
expressions (5.2) and (5.5) for the general law.

5.2.1 Mass conservation

For a system (material form):

Mass conservation requires the mass of any system (material domain) to be constant, since
there is no material flow through the boundary of the system (mass conversion to energy is
not considered here). Therefore, the material time derivative of mass must vanish, i.e.:
D
Dt 
dV = 0 . (5.8)
V
The rate of change of
the mass in the system

The differential form of the mass conservation law or the continuity equation can be obtained
using the following form of expression (5.8):

 (,t0 )dV0 =  ( x,t)dV , (5.9)

V0 V

where V is the volume occupied by a material at time t that was occupying volume V0 at t0,
and 0 =  ( , t0) and  =  (x, t) are the densities at t0 and t, respectively. By substituting
variables (x = x ( , t)) in the integral on right hand side of (5.9), one gets:

 ( x,t)dV =  ( x(,t),t )JdV0 , (5.10)

V V0

where J is the Jacobian of the transformation x = x ( , t). It follows that:

 (0  J )dV0 = 0 , (5.11)

V0

54
 5. FUNDAMENTAL LAWS

and since V0 is arbitrary volume:

J =  0 , (5.12)
or:
D D
Dt
( J ) = 0 or Dt
(dV ) = 0, (5.13)

since dV = J dV0 and dV0 = const. This equation is Lagrangian or material form of the
continuity equation. From (5.13) it follows that:
D D D D
Dt
( dV ) =  ( dV ) + 
Dt Dt
dV = 
Dt
dV , (5.14)

=0
which has an important implication for other conservation laws.

For a control volume (spatial form):

The mass conservation law for a control volume takes the form:


 t dV +  v ndS = 0 , (5.15)
V S
Change of mass Mass flux
in the CV through CS

Based on (5.7), the differential continuity equation in the spatial form can be written in the
form:

+ div( v) = 0 . (5.16)
t
Local Convective
change change

or:
 D
+ v grad  + div v = + div v = 0. (5.17)
t Dt

According to (2.76), (4.62) and (5.16), it follows that:

D(  )  (  )  
= + div (v) =  +  + v  grad  +  div(v ) =
Dt t t t
, (5.18)
      D
 + div(v ) +   + v grad  = 
t t Dt
D
=0 =
Dt
which means that D()/Dt = D/Dt can be expressed in two different forms:

55
 5. FUNDAMENTAL LAWS

D(  ) D (  )   
= = + div (v) =   + v grad  . (5.19)
Dt Dt t t
Strong conservation Weak conservation
form form

5.2.2 Reynold’s theorem for a density-weighted integrand

The material time derivative of an integral (4.76) or (4.77) can be written in the form:
D D
Dt  FdV = 
Dt
( FdV ) , (5.20)
V V

since by using (4.74) it can be shown that:

D DF D(dV ) DF
( FdV ) = dV + F = dV + Fdiv v dV . (5.21)
Dt Dt Dt Dt
If integrand F in (5.20) is replaced by F and recalling (5.14), expression (5.20) can be
written as:
D D DF
 FdV =  (FdV ) =   Dt dV . (5.22)
Dt Dt
V V V

This is a consequence of Reynold’s transport theorem and mass conservation, and is useful
expression in continuum mechanics. For example, one can see that (5.8) and (5.13) follow
directly from (5.22).

5.2.3 Conservation of linear momentum (1st Euler’s law, 2nd Newton’s law)

For a system (material form):

The (linear) momentum conservation principle states that the material time derivative of the
linear momentum is equal to the resultant force acting on the system (also see (5.2)), i.e.:
D
Dt   fbdV   ndS .
vdV = + (5.23)
V V S
Rate of change of Total body Total surface
total linear momentum force force

According to (5.4) and (5.18), the differential form of the momentum equation takes the form:
Dv
 = fb + div  . (5.24)
Dt
Rate of change of Body Surface
linear momentum force force

For a control volume (spatial form):

The momentum conservation law for a control volume takes the form:

56
 5. FUNDAMENTAL LAWS

 (v )
 t
dV +  v( v n)dS =  fbdV +    ndS , (5.25)
V S V S
Change of mom. Flux of mom. Total body Total surface
in CV across CS force force

and according to (5.7), its differential form is:

( v)
+ div (vv ) = fb + div  . (5.26)
t
Local Convective. Body Surface
change change force force

Equations (5.24) and (5.26) can be expressed in the index form in Cartesian coordinate
system as:
Dvi  ji
 = fb,i + , (5.27)
Dt x j

which can be written in the expanded form as:

Dv x   yx  zx
 = f b,x + xx + + ,
Dt x y z
Dv y  xy  yy  zy
 = f b,y + + + , (5.28)
Dt x y z
Dv   yz  zz
 z = f b,z + xz + + .
Dt x y z

Equations (5.27) or (5.28) are known as Cauchy equations.

5.2.4 Conservation of angular momentum (2nd Euler’s law)

For a system (material form):

The integral form of the conservation of angular momentum is obtained by taking the cross-
product of each term in the corresponding linear momentum principle with the position
vector x. Hence, it is also called the moment of momentum principle, and is expressed as:
D
Dt 
 (x  v)dV =  ( x  fb )dV +  ( x   ) ndS . (5.29)
V V S
Rate of change of tot. Total moment Total moment
moment of momentum of body forces of surface forces

According to (5.4), the differential form of the moment of momentum principle takes the
form:

57
 5. FUNDAMENTAL LAWS

D( x  v)
 = x  fb + div (x   ) . (5.30)
Dt
Rate of change of Moment of Moment of
moment of momentum body forces surface forces

This is also known as the second Cauchy law. The first term on the left and the second term
on the right hand side of (5.30) can be expressed in the following way:
D( x  v) Dv Dx Dv Dv
=x +v =x +v v= x  ,
Dt Dt Dt Dt Dt (5.31)
div( x   ) = x  div  + grad x  .

Based on (5.31), expression (5.30) can be written as:

 Dv 
x    fb  div  = grad x   . (5.32)
Dt

and according to (5.24) and (2.77) it follows that:

I   = 0 or ijk kj = 0 , (5.33)

hence:
 T =  or  ji =  ij . (5.34)

This means that the second Cauchy law can also be formulated in the following way: the
necessary and sufficient condition for the moment of momentum principle to be satisfied is
the symmetry of the Cauchy stress tensor.

For a control volume (spatial form):

The moment of momentum principle for a control volume takes the form:

[ ( x  v)]
 t
dV +  (x  v)( v n)dS =  ( x  fb )dV +  (x   ) ndS . (5.35)
V S V S
Change of mom. of Flux of mom. of Total moment of Total moment of
momentum in CV momentum across CS body forces surface forces

and according to (5.7), its differential form is:

[ (x  v)]
+ div[ v( x  v)] = x  fb + div (x   ) . (5.36)
t
Local Convective. Moment of Moment of
change change body forces surface forces

58
 5. FUNDAMENTAL LAWS

5.2.5 Energy conservation law (1st law of thermodynamics)

For a system (material form):

Here, only thermodynamical processes where the sources of energy are mechanical work and
heat are considered. The principle of conservation of energy, i.e. the energy balance
principle, states that the rate of change of total energy is equal to the work done by surface
and body forces plus the heat delivered to the body by the heat flux and other heat sources.
The integral form of the energy balance principle takes the form:
D
Dt 
edV =   q ndS +  hdV +  fb  vdV +  (  v)  ndS , (5.37)
V S V V S
Rate of change of Heat Heat Work done (power)
total energy delivered source

where the sign of the heat flux term is negative since positive heat flow is out of the body.
Here, q is the heat flux, h is the heat source, and e is the total specific energy. In continuum
mechanics, e is often taken as the sum of internal energy u and kinetic energy, i.e.:
1 v2
e = u + vv= u + . (5.38)
2 2

According to (5.4) and (5.18), the differential form of the energy balance equation takes the
form:
De
 =  divq + h + fb  v (  v) .
+ div (5.39)
Dt
Rate of change Heat Heat Work by Work by
of energy flux source body forces surface forces

For a control volume (spatial form):

The energy conservation law for a control volume can be expressed as:
( e )
 dV +  e( v n)dS =  q ndS +  hdV +  fb  vdV +  (  v)  ndS ,(5.40)
t
V S S V V S
Change of Energy flux Heat Heat Work done (power)
energy in CV through CS delivered source

and according to (5.7), its differential form is:

 De
 =  Dt

(5.41)
(e)
+ div (ev) = divq + h + fb  v + (  v)
div
t
Local Convective Heat Heat Work by Work by
change change flux source body forces surface forces

59
 5. FUNDAMENTAL LAWS

Equation of mechanical energy

A dot product of the momentum equation and the velocity vector:

Dv

Dt
v = ( fb + div )  v (5.42)

gives the so called equation of mechanical energy:

D v 
2
 = fb  v + (div )  v .
Dt  2

(5.43)

Change of Work by Mechanical

Mech. energy body forces work

It is important to notice that this equation is a consequence of the momentum equation and
not the energy equation.

Heat equation
Using (5.38), the energy equation (5.39) [or (5.41)] can be written in the form:
Du Dv
 + v  ( )  v + :gradv ,
= divq + h + f b  v + div (5.44)
Dt Dt

giving:
Du  Dv 
 + v    fb  div() = divq + h + :gradv . (5.45)
Dt Dt

By employing (5.24) [or (5.26)], one gets:

Du
 =  divq + h +  :gradv . (5.46)
Dt
Total change Heat Heat Deformation
of internal energy flux source work

This equation is known as the heat equation. By summing the equation of mechanical energy
and the heat equation one gets the equation of total energy (5.39) [or (5.41)]. Equation (5.46)
in the index notation takes the form:
Du q j v j
 = + h +  ij , (5.47)
Dt x j xi

and its expanded form is:

60
 5. FUNDAMENTAL LAWS

Du  q qy qz  v x v v
 = 
x + +  + h +  xx +  yx x +  zx x +
Dt  x y z x y z
v y v y v y
 xy +  yy +  zy + (5.48)
x y z
v v v
 xz z +  yz z +  zz z
x y z

In deriving (5.46), the following expression was used [see (2.76)]:

(  v)
div = ( )v
div +  :gradv . (5.49)
Total work by Motoric Deformation
surface forces work work

While the motoric work affects the mechanical energy change, the deformation work is
converted into heat and increases the internal energy in the control volume. By decomposing
the stress tensor into its spherical and deviatoric parts (3.26), one gets:

( d
:gradv =  pI +  :gradv = )  pdivv +  d :gradv. (5.50)
Deformation Reversible Irreversible
work part part

It can be seen that the deformation work consists of the reversible part (which is a
consequence of the spherical stress component) which can be converted back into the
mechanical energy, and of the irreversible part (which is a consequence of the deviatoric
stress component) which cannot be converted back into the mechanical energy.
When the heat flux and heat sources vanish, i.e. in a purely mechanical process, the energy
equation becomes [see (4.81) and (4.82)]:
Du
 =  :gradv =  :(D  V) =  :D =  ij D ji , (5.51)
Dt
which is no longer a partial differential equation. It can be seen that the internal power (the
internal energy rate) is given by the contraction (scalar, inner or double-dot product) of the
rate of deformation and the Cauchy stress. We therefore say that the rate of deformation and
the Cauchy stress are conjugate in power (or in work or in energy).

5.2.6 Law of entropy production (2nd law of thermodynamics)

The first law of thermodynamics does not answer the question of the extent to which the
energy conversion process is reversible or irreversible. The second law of thermodynamics
gives the basic criterion for irreversibility with its statement on the limitations of entropy
production. In continuum mechanics the specific entropy (per unit mass) or entropy density is
usually denoted by s, so that the total entropy is  sdV . The entropy can be viewed as the
V

61
 5. FUNDAMENTAL LAWS

quantitative measure of microscopic randomness and disorder. The entropy of a system can
change either through interactions with surroundings, or by changes that take place within the
system, i.e.:
ds = dse + dsi , (5.52)

where ds is the increase in specific entropy, dse is the increase due to interaction with the
exterior, and dsi is the internal increase. The change dsi is never negative, it is zero for a
reversible process, and positive for an irreversible process:
dsi > 0 (irreversible process)
. (5.53)
dsi = 0 (reversible process)

In a reversible process, if dqR denotes the heat supplied per unit mass to the system, the
change dse is given as:
dq R
dse = (reversible process), (5.54)
T
where T is the absolute temperature.

For a system (material form):

According to the second law, the rate of increase of total entropy of the system is never less
than the sum of the entropy influx through the surface plus the entropy produced internally by
body sources. Mathematically, this entropy principle may be expressed in integral from as the
Clausius-Duhem inequality:
D 1
Dt  sdV   q ndS +
T  dV , (5.55)
V S V

where  is the local entropy source per unit mass.  is sometimes given as  = r / T , where
r is internal heat supply per unit mass and time (e.g. from radiation). The equality in (5.55)
holds for reversible processes, while the inequality applies to irreversible processes. The
difference between the rate of change of entropy and the rate of entropy input in a body
determines the total production of entropy , i.e.:
D 1
Dt 
= sdV +  q  ndS   dV  0. (5.56)
T
V S V

This relation describes the direction of the energy transfer and postulates irreversibility of
various thermodynamic processes (observations show that heat flows from the warmer to the
colder regions of a body; mechanical energy can be transformed into heat by friction, and this
can never be converted back into mechanical energy). The situation  < 0 never occurs; it
would mean that molecules organise themselves globally. A thermodynamic process is called
reversible if it is not accompanied by any entropy production ( = 0). Real processes are

62
 5. FUNDAMENTAL LAWS

always irreversible and characterised by  > 0. These are always associated with energy
dissipation.

For a control volume (spatial form):

The Clausius-Duhem inequality for a control volume can be expressed as:

( s) 1
 dV +  s(v  n)dS    q  ndS +  dV . (5.57)
t T
V S S V

The differential form of the expression (5.55) [or (5.57)] is:

 Ds 
 =  Dt
(5.58)
(s) q
+ div (sv)   div + 

t T
By using (2.76):
q 1 1
div = divq  2 q gradT , (5.59)
T T T
giving:
Ds 1 1
   divq + 2 q gradT +  . (5.60)
Dt T T
1
Since the heat flows from hot to cold, we have  2 q gradT  0 . Expression (5.58) gives
T
the local form of the entropy production rate  per unit mass as:
Ds 1 q
 + div    0. (5.61)
Dt  T

This inequality must be satisfied for any process, and for this reason it plays an important role
in imposing restrictions upon constitutive relations discussed later.
From Gibbs relation (for an ideal gas) we have that:
 1
du = Tds  pdv = Tds  pd , (5.62)
 

where v is specific volume, and this can be written in the following form:
Ds Du D(1/ )
T = p . (5.63)
Dt Dt Dt
Using the continuity equation (5.17), the last term in the above equation becomes:
D(1 /  ) 1 D 1
= 2 = divv. (5.64)
Dt  Dt 

Now, with the help of (5.64) and (5.46), expression (5.63) takes the form:

63
 5. FUNDAMENTAL LAWS

Ds
T =  divq + h +  :gradv + pdivv, (5.65)
Dt
and by using (5.50) and the fact that  p =  p (3.25), one gets:
Ds
T =  divq + h  pdivv +  d :gradv + pdivv =  divq + h +  d :gradv . (5.66)
Dt
From (5.59) it follows that:
1 q 1
divq = div + 2 q  gradT , (5.67)
T T T
and finally, (5.66) takes the form:
Ds q 1 h 1 d
 =  div  q gradT + +  :gradv . (5.68)
Dt T T2 T T
Total increase Reversible Increase in entropy Increase in entropy
of entropy heat exchange because of irreversible because of irreversible
heat exchange due to exchange of mechanical
gradT work to heat

For an adiabatic process with divq + h = 0 , we derive that:

Ds 1 d
 =  :gradv .
Dt T

Expression (5.68) written in the index notation becomes:

Ds  q j  1 T h 1 d v j
 = 
  2 qj + +  , (5.69)
Dt x j  T T x j T T ij xi

and its expanded form is:

Ds  
q x  
qy  
qz  
 =  +   +   
Dt  x  T y  T z  T 
1
T T T  h
2  q x x + qy y + qz z  T +
+
T
1 v v x v 

T
(  xx + p) x +  yx
x y
+  zx x +  .
z
(5.70)

 v y v y v y 

(
 xy x +  yy + p ) y
+  zy +
z 
 v z v z v z 
 xz x +  yz y + ( zz + p) z 


For an irreversible adiabatic process Ds/Dt>0 by the second law, and in much of continuum
mechanics, it can be shown from relations derived similarly to (5.68) that the dissipation
d
function  :gradv must be positive, since both  and T are positive.

64
 5. FUNDAMENTAL LAWS

5.2.7 Discussion
Five fundamental laws of continuum mechanics presented above, result in a system of 9
balance equations in Cartesian coordinate system: 1 continuity, 3 momentum, 1 energy, 3
moment of momentum, and 1 entropy equation. This system, however, is an open system as it
contains 19 unknown scalar functions: , vj, ij, qi, u, T and s (i = 1, 2, 3 and j = 1, 2, 3). To
close the system, 10 additional equations are required, which will link the 19 dependent
variables. For this purpose, the following constitutive equations are used:

1) Relationship between stress tensor and deformation tensor and/or rate of deformation
tensor (6 equations),
2) Relationship between the heat flux and temperature gradient (3 equations),
3) State equation (1 equation).
Basic constitutive equations will be given in the next chapter.

5.3 Equations for large deformations

We have seen that there are two basic descriptions of motion and deformation: Lagrangian
and Eulerian. Similarly, one can define different stress measures with respect to the reference
and current configurations. The Cauchy stress tensor, defined in Chapter 3 as a function of
spatial coordinates x in a current configuration, is the only physically meaningful stress
quantity. When the deformation is defined in terms of material coordinates  in a reference
configuration, then the stress should also be expressed as a function of material coordinates,
and the balance equations defined in a reference configuration. These are often called
Lagrangian conservation equations.
Lagrangian formulation is convenient for elastic solid mechanics problems in particular,
because the elastic body is assumed to return back to its reference undeformed state after the
removal of forces (mechanical and thermal). However, the equations of motion (or
equilibrium equations) must be satisfied in the deformed or current configuration, and the
corresponding stress is also defined in the deformed configuration (Cauchy stress). Therefore,
both the stress and the deformation should be defined in the same configuration either
deformed or undeformed.

An important condition when choosing mechanical quantities for a constitutive relationship is

that they are complementary. Specifically, the stress and deformation measures must be
conjugate. This means that an increment of deformation multiplied (in a certain way) with the
stress must be equal to an incremental work per unit volume of material [see section 5.2.5,
eq. (5.51)]. Therefore, material’s volume, surface and length must be consistent with the
chosen stresses and deformations.

65
 5. FUNDAMENTAL LAWS

5.3.1 Kirchoff and Piola-Kirchoff stress tensors

Apart from Cauchy stress tensor  , defined using Eulerian description (Chapter 3), which
relates the current force with the current surface:
df = n  dS , (5.71)

there are three additional stress tensors, also called pseudo-stress tenors, defined as:
1. Kirchoff stress tensor  (1) obtained by scaling the Cauchy tensor:
(1) 
 = 0  = J . (5.72)


2. First Piola-Kirchoff stress tensor  (1), gives the current force acting on the deformed
configuration per unit of the undeformed area dS0:
(1)
df = n0  dS0 . (5.73)

3. Second Piola-Kirchoff stress tensor  (2), gives the current force transformed to the
reference configuration df0 per unit of the undeformed area, where df = F df0 :

df0 = n0  (2)dS0 , or df = F n0  (2)dS0 , (5.74)

where the relationship between the deformed area dS = ndS and undeformed area
dS0 = n0dS0 is given by Nanson’s relation:

ndS = 0 n0 F 1dS0 = Jn0 F 1dS0 . (5.75)


Forces and surfaces involved in defining Piola-Kirchoff tensors are illustrated in Fig. 5.1.

df 0

n0 df

dS 0
df
dS
n
S V
S0 V0

Figure 5.1 Surface forces and surface elements involved in defining Piola-Kirchoff tensors

Using expressions (5.71), to (5.75), one can relate Cauchy with Piola-Kirchoff stress tensors.
First, the current force acting on the current surface can be expressed as:

66
 5. FUNDAMENTAL LAWS

df = n  dS = n0  (1)dS0 = F n0  (2)dS0 , (5.76)

and it follows that:

 (1)  
=  = F   (1) = F   (2)  FT , (5.77)
0 0 0

(1)  1 (2) T
 = 0 F  =   F , (5.78)



( ) ( )
T T
 (2) = 0 F 1    F 1 =  (1)  F 1 . (5.79)


These expressions show that second Piola-Kirchoff stress tensor is symmetric, while the first
Piola-Kirchoff stress tensor is not, and is therefore inconvenient to use.

5.3.2 Momentum equation

In this section we derive the momentum equation from the condition that the rate of change
of momentum of a material in V, which occupied V0 in the reference (undeformed)
configuration at t0, is equal to the sum of all forces acting on the material in V. It is important
to notice that in what follows, although the integrals are calculated over S0 and V0, the forces
are current forces acting on the material in V and through S, and the momentum is current
momentum of the material in V at time t.

Integral form:
If substitution x = x(,t) is introduced in the equation of motion (5.23), i.e. the integrals are
expressed in the reference configuration, then:
D
Dt  vJdV0 =  fb0 JdV0 +  n0  (1)dS0 . (5.80)
V0 V0 S0

Taking into account that  0 = J (5.12), one gets the momentum equation for the reference
configuration in terms of the first Piola-Kirchoff stress:
D
 0vdV0 =  0fb0dV0 +  n0 
(1)
dS0 , (5.81)
Dt
V0 V0 S0

where fb0 (,t) = fb[ x(,t)] . This expression is completely analogous to (5.23) written for the
current configuration with Cauchy stress tensor. Since the first Piola-Kirchoff stress tensor is
not symmetric, the second Piola-Kirchoff stress is more often used in the momentum
equations for the reference configuration:
D
Dt  0vdV0 =  0fb0dV0 +  n0  (2)  FT dS0 , (5.82)
V0 V0 S0

67
 5. FUNDAMENTAL LAWS

which is obtained by using (5.81) and (5.78).

Differential form:
The differential form of equations (5.81) and (5.82) is:
Dv (1)
0 =  0fb0 + div  , (5.83)
Dt

0
Dv
Dt (
=  0fb0 + div  (2) F T , ) (5.84)

while the corresponding equations in the index form in Cartesian coordinate system are:

Dvi  (1)
ji
0 = 0 fb 0,i + , (5.85)
Dt  j

Dvi   (2) x i 
0 = 0 fb 0,i + 
 j
 jk k 
. (5.86)
Dt

5.3.3 Moment of momentum equation

The moment of momentum equation in the reference configuration reduces to the symmetry
of the second Piola-Kirchoff stress  (2).

5.3.4 Energy equation

The deformation work per unit time (power) V :gradv dV = V :D dV [see (5.46) and
(5.51)] in volume V of the current configuration can be written for the volume V0 occupied by
the same material in the reference configuration using the first Piola-Kirchoff stress tensor as:

Pdef =  :gradv dV =  F  (1) :gradv JdV0 =  F   (1) :gradv dV0 =
0
V V0 V0
. (5.87)
x (1) v j (1) v j x
 ik kj xi dV0 =  kj x i ik dV0
V0 V0

Since:
DF D  x j   Dx j v j v j xi
= 
i ji k = 
i ji k = i j ik = i i = Y F, (5.88)
Dt Dt  k k  Dt k xi  k j k

then, expression (5.87) can be written in the form:

T
(1) DF
Pdef =   :
Dt
dV0 . (5.89)
V0

68
 5. FUNDAMENTAL LAWS

If the deformation work is expressed by V :D dV and the second Piola-Kirchoff stress

tensor is used, then with the help of (5.77) and (5.12) one gets:

Pdef =   :D dV = 

0( ) (
F  (2)  FT :D JdV0 =   (2): F T D  F dV0 , ) (5.90)
V V0 V0

DL T
and since = F  D F , where L is defined by (4.30), we have that:
Dt
DL 1 (2) DG
Pdef =   (2) : dV0 =   : dV0 . (5.91)
Dt 2 Dt
V0 V0

The total heat delivered through the current area S can be written in terms of the undeformed
surface S0 by using Nanson’s relation (5.75):

(
  q  ndS =   q  Jn0 F 1dS0 =   JF 1 q  n0dS0 , ) (5.92)
S S0 S0

where the term in brackets can be regarded as the heat flux per unit reference area. Finally,
the heat equation in the integral form can be written for the reference configuration with the
first Piola-Kirchoff stress tensor as:
DFT
D
Dt  ( )
0 udV0 =   JF 1  q  n0dS0 +  h0dV0 +   (1) : Dt
dV0 , (5.93)
V0 S0 V0 V0

or with the second Piola-Kirchoff stress tensor as:

D
Dt  0 udV0 ( )
=   JF 1  q  n0dS0 +  h0dV0 +   (2): Dt
DL
dV0 , (5.94)
V0 S0 V0 V0

where h0 = Jh is the heat source per unit initial volume.

5.4 Variational methods

Variational methods play an important role in continuum mechanics. The fundamental
variational principle is the principle of virtual displacement or principal of virtual work. This
is a fictitious work calculated using statically admissible forces. These forces are assumed to
remain unchanged while they do work on infinitesimal kinematically admissible
displacements. Moreover, the forces and displacements can be prescribed independently.

5.4.1 Definition of a mechanical problem

The principle of virtual work is an alternative way of expressing equilibrium equations in
classical mechanics. It will be shown that this equivalency holds in continuum mechanics as
well.

69
 5. FUNDAMENTAL LAWS

Boundary conditions
Let’s consider a mechanical problem in region V with boundary conditions prescribed on S.
The problem is said to be well posed if at every point on S and for every direction of a
coordinate system (e.g. local coordinate system with one axis normal to the surface S), there
is one and only one condition prescribed via either displacement components ui (or velocity
vi) or surface tractions fsi:
ui = ubi or f si = fbi (i = 1,2,3) on S . (5.95)

Virtual displacements

The virtual displacement field u is said to be kinematically admissible if it satisfies any

prescribed displacement boundary conditions. It must also be continuous and it must possess
continuous first partial derivatives in V. Since the virtual displacements are additional
displacements from the equilibrium configuration, a corresponding virtual displacement
component must be zero in every point of S where the actual displacement component is
prescribed, i.e.:
ui = 0 if ui = ubi (i = 1,2,3) on S . (5.96)

s
Virtual strain tensor E corresponds to the virtual displacement u , and its components are
given as:

 ij =  +
( )
.
1   (ui )  u j
(5.97)
2  x j x i

Virtual forces
When certain external (body and surface) forces are acting on a body, a stress state is said to
be statically admissible if stress tensor 

1. Satisfies the equilibrium equations [eqs. (5.24) or (5.26) with Dv/Dt = 0] in the interior V
of the body:

 fbdV +   ndS = 0 . (5.98)

V S

2. Satisfies boundary conditions on S wherever the boundary tractions are prescribed:

n j ji = fbi if fsi = f bi (i = 1,2,3) on S . (5.99)

70
 5. FUNDAMENTAL LAWS

5.4.2 Virtual work principle

Principle

If a problem is well posed, virtual displacements kinematically admissible and virtual forces
statically admissible, then the virtual work of external forces is equal to the virtual work of
internal forces:

 fb  u dV +  n   u dS =   :Es dV . (5.100)

V S V
Virtual work Virtual work
of external forces of internal forces

This condition is equivalent to the equilibrium equations.

Proof of equivalency
Let us show that the virtual work principle (5.100) follows from the equilibrium equation
(5.98). Since forces remain unchanged during virtual displacement, then the dot product of
the total external force and virtual displacement u , gives the virtual work of external forces:
Wext =  fb  u dV +   n  u dS . (5.101)
V S

By applying the Gauss theorem to the surface integral in (5.101), one gets:
Wext =  [:gradu + ( fb + div)  u ]dV =  :E s dS , (5.102)
V S

because the term in the curly bracket is zero (equilibrium equation), symmetric part of
:grad u is equal to the virtual strain tensor Es , and the double dot product of symmetric
stress tensor  and anti-symmetric part of the tensor :grad u is zero.

Expressed in terms of Piola-Kirchoff tensors

It can be shown that the virtual work principle (5.100) in reference configuration expressed in
terms of the first Piola-Kirchoff stress tensor takes the form:

 0fb 0  x dV0 +  n0  (1)  x dS0 =   (1) :FT dV0 . (5.103)

V0 S0 V0

where u = x  , so u = x , and in terms of the second Piola-Kirchoff stress tensor:

 0fb 0  x dV0 +  n0  (2) F T  x dS0 =   (2) :L dV0 . (5.104)

V0 S0 V0

(1) T (2)
From the above two expressions, one can see that  and F , and  and L are work
conjugate pairs.

71
 5. FUNDAMENTAL LAWS

Generalisation

It can also be shown that the virtual work principle (5.100) is valid if the inertial force is
added to the body force. Namely, by taking that fb + div = Dv / Dt in (5.104), one gets:
 Dv 
 fb  u dV +  n   u dS =   :E
s
+  u
dV , (5.105)
Dt
V S V

or
1
 fb  u dV +  n   u dS =   :E dV +    v 2 dV .
s
(5.106)
2
V S V V
Work by external forces Deformation Change of kinetic
Work energy

In deriving (5.106), the following was used:

Dv Dv 1 Dv 2
 Dt
 u dV =  
Dt
 vt dV =  
2 Dt
t dV
V V V
. (5.107)
 1  1
=    v 2 dV  t =    v 2 dV
t  2 

2
V V

Expression (5.106) represents the equation of mechanical energy:

K = Wext  Wdef , (5.108)

where K is the total kinetic energy, and Wdef is the total work of deformation. By inserting
(5.108) into the energy balance (first law of thermodynamics)

U + K = Q + Wext , (5.109)
=E

one gets the heat equation:

U = Q + W def , (5.110)

where E is the total energy and U is the total internal energy. Therefore, the virtual work
principle can be considered as an energy principle.

72
 6. CONSTITUTIVE RELATIONS

Chapter 6

Constitutive relationships
6.1 Introduction

6.1.1 Ideal materials

So far, we have considered quantities, which characterise the motion and stresses of any
continuum, and we have derived the links between them (fundamental laws) independently of
the physical properties of a continuum. However, different materials behave differently under
loading (e.g. rubber specimen will deform differently under the same load than the steel
specimen of the same dimensions). Unlike solids, fluids and gasses will change their shape
under relatively small forces. Gasses, unlike fluids, will considerably change their volume
under small changes in pressure. The physical properties of a continuum must be taken into
account when analysing the motion of that continuum.
In this chapter, we shell consider mathematical relations between static (stresses), kinematic
(displacements, deformations, velocities) and thermal (heat flux, temperature) variables, which
define the behaviour of a continuum subjected to mechanical and/or thermal loading. These
relations are called constitutive relations, as they describe macroscopic behaviour resulting
from the internal constitution (structure) of the material. But materials, particularly in the solid
state, behave in extremely complex ways when a whole range of possible temperatures and
deformations are considered. Therefore, it is not feasible (neither it is possible at this stage) to
describe accurately a real material over its entire range of behaviours. Instead, we formulate
separately constitutive equations describing various kinds of ideal material response, each
designed to approximate material’s response over a suitably restricted range.

There is a large number of different materials, hence there is a large number of corresponding
constitutive equations. However, it is somewhat surprising that two constitutive relationships,
one for linear elastic solid and another for linear viscous fluid, can well describe behaviour of
many real materials for a range mechanical and thermal loading. Of course, real materials are
always more or less different than the two simple ideal materials, and when the differences

73
 6. CONSTITUTIVE RELATIONS

become considerable, then one is dealing with non-Newtonian fluids, viscoelastic, plastic or
viscoplastic solids, and corresponding constitutive relations are much more complex.
The ideal elastic (Hookean) solid has the property that it returns to its initial (undeformed)
state after the removal of all mechanical and/or thermal loads. Its constitutive equation is given
by a linear relation between the stress tensor and strain (deformation) tensor of the form:
Stress tensor = f(Deformation tensor).

The ideal Newtonian fluid differs from a solid in that it cannot support a shear stress; in a fluid,
shear deformation (flow) will continue as long as any shear stress is applied. Its constitutive
equations is expressed by a linear relation between the stress and the rate of deformation
(strain rate) tensor of the form:
Stress tensor = f(Rate of deformation tensor).

Unlike the elastic material, viscous fluid has no memory of its initial state; the stress at the
point depends exclusively on the instantaneous rate of deformation at that point.

Not all material can be separated into these two categories (where the stress is only a function
of either deformation or the rate of deformation). The elastic solid and viscous fluid are the
two extreme cases of a wide range of materials including, viscoelastic solids, elastic solids that
have rate dependence but still possess well defined initial undeformed state to which they
eventually return when the loads are removed, or viscoplastic materials, which remind more of
a thick fluid but are still exhibiting some aspects of a solid, etc.

6.1.1 Isotropic tensors

It was said in Section 2.3 that a material is isotropic with respect to a certain property, if that
property is the same in all directions. Physical properties of materials are usually described
using tensors (of various order), and the precise definition of material isotropy can be given by
defining an isotropic tensor.
A tensor is isotropic if its Cartesian components are unchanged by any orthogonal
transformation of the coordinates. A trivial example is the zero tensor of any order. All tensors
of zero order (scalars) are isotropic, but there are no isotropic first-order tensors (vectors)
except the zero vector. The unit tensor I, whose components are given in any Cartesian system
by Kronecker delta ij, is isotropic, and it can be shown that I and sI, where s is a scalar, are
the only non-trivial isotropic second-order tensors. The only non-trivial isotropic third-order
tensors are the ones whose Cartesian components are given by the permutation symbol ijk, and
scalar multiples of it.
It can be shown that the Cartesian components of all isotropic fourth-order tensors can be
written in the form:
( ) ( )
cijkl =  ij kl + μ ik  jl + il  jk +  ik  jl  il jk , (6.1)

74
 6. CONSTITUTIVE RELATIONS

where , μ and  are constants which have the same value in all coordinate systems. If cijkl is
symmetric in either ij or kl, the last bracketed term in (6.1) is zero and we have:
( )
cijkl =  ij kl + μ ik  jl + il  jk , (6.2)

which means that any isotropic fourth-order tensor cijkl symmetric in ij or kl is described with
only two constants (scalars)  and μ.

6.2 Classical constitutive relations and equations of state

Classical constitutive relationships are related to situations with very simple states of stress and
deformations or flow. They are:

1. Hooke’s law of elasticity,

2. Pascal’s law of hydrostatic pressure,

3. Newton’s law of viscosity,
4. Fourier’s law of heat conduction.

These laws describe the simplest, so called linear materials. They are later generalised to give
constitutive equations for arbitrary states of stress and deformation (or flow), and are able to
describe reasonably well the behaviour of a large number of real materials.

6.2.1 Hooke’s law

By analysing results of his tensile tests on metals, Hooke established that, away from the ends
where the specimen was loaded (Fig. 6.1), the following relation holds:
F l
=E , (6.3)
S l
or:
 = E , (6.4)

where
is the normal stress,  is the strain (relative extension), and a proportionality constant
E is so called Young’s modulus.

F S
F

Figure 6.1 Tensile test

75
 6. CONSTITUTIVE RELATIONS

6.2.2 Pascal’s law

Pascal experimentally established that a stationary fluid acts on a stationary solid surface by a
force which is always normal to the surface, and that its intensity does not depend on the
orientation of the surface. The ratio between the force intensity and the surface:
F
p=  fs =  pSn , (6.5)
S
was named hydrostatic pressure.

6.2.3 Newton’s law of viscosity

Studying the fluid flow between two parallel close surfaces (Fig. 6.2), Newton came to the
conclusion that the force required to maintain the flow is proportional to the surface and the
gradient of the velocity, i.e.
F v
=μ , (6.6)
S h
or expressed in the differential form:
dv
 =μ , (6.7)
dy

where  is the tangential (shear) stress, dv/dy is the velocity gradient, and the proportionality
coefficient μ is the coefficient of (dynamic) viscosity.
v F

h
y

Figure 6.2 Fluid flow between the two parallel planes

6.2.4 Linear and non-linear materials

The common feature of the Hooke’s law, which is applicable for solids, and the Newton’s law,
which is applicable for fluids, is that the relationship between the stresses and strains or the
rate of deformation is linear. Therefore, materials whose behaviour can be described by these
laws are called linear materials.

Most solids subjected to relatively small stresses and low temperatures behave according to the
Hooke’s law. However, this is no longer the case under substantial loads and/or high

76
 6. CONSTITUTIVE RELATIONS

temperatures, where materials behave in inelastic manner. This inelasticity can be manifested
via plastic flow, hardening, creep, relaxation, etc.
Similarly, almost all gasses and low-density fluids behave according to the Newton’s law.
However, there are some important fluids with non-linear relation between the stresses and rate
of deformation. These are called non-Newtonian fluids.

Inelastic solids and non-Newtonian fluids are also classified as non-linear materials.

6.2.5 Fourier’s law of heat conduction

Fourier established that the amount of heat per unit time and per unit area conducted between
two surfaces, which are at constant temperatures T1 and T2 and distance h apart, is:
T T
q = k 2 1, (6.8)
h

or in the differential form:

dT
q=k , (6.9)
dx
where k is the coefficient of thermal conductivity. It is important to notice that this law is
applicable for all materials, unlike the previous ones which are specific for solids or fluids.

6.2.6 Caloric equation of state

In the state of thermodynamic equilibrium, the (local) specific internal energy u in a given
point of a system depends on the temperature T, and some (not necessarily all) aspects of
deformation of the particle from its reference state to the current state which coincide with the
point. So, u is said to be determined by the thermodynamic state, specified by a number of
substate variables and the specific entropy s. The number of substate variables depends on the
type of the material. For example, in the simplest case of a fluid pure substance (an ideal
viscous fluid) there is only one substate variable, specific volume v = 1/  [see (5.62)]. In ideal
elasticity, there are 6 substate variables, the components of one of the strain or deformation
tensors.

For a viscous fluid we have that:

 u 
u = u(T, v )  du =   dT = c v dT , (6.10)
T v

where cv is the specific heat at constant volume, while for an elastic solid we have:
 u
(
u = u T, ij )  du =  dT = c E dT ,
T 
(6.11)
ij

77
 6. CONSTITUTIVE RELATIONS

where cE is the specific heat at constant deformation. In the theory of linear elasticity, either cv
and cE can be used in (6.11), as their values are close to each other.

6.3 Elastic solids

A material is called ideally elastic when a body recovers its original form upon removal of all
forces, and if there is a one-to-one relationship between the state of stress and the state of
strain. If the components of the displacement gradient are small compared to unity, then the
differences between the current (deformed) and initial (undeformed) configurations can be
neglected, i.e. no distinction need be made between material ( ) and spatial (x) coordinates.

Under the above conditions it can be assumed that the Cauchy stress tensor is proportional to
the small strain tensor:
 = C:E s or  ij = cijkl kl , (6.12)

where C is the fourth-order tensor with 81 components. Since  and Es are symmetric tensors,
then:
cijkl = c jikl and cijkl = c ijlk , (6.13)

and number of (independent) components reduces to 36.

6.3.1 Generalised Hooke’s law

If a material is elastically isotropic, i.e. the elastic constants are the same for all possible
choices of Cartesian coordinates, then cijkl are the components of the fourth-order isotropic
tensor symmetric in the indices ij and kl. Based on (6.2), expression (6.12) can now be written
in the form:

[ ( )]
 ij = ij  kl + μ ik jl + il  jk kl , (6.14)

whence we obtain generalised Hooke’s law:

( )
 ij =  kkij + 2μ ij or  =  trE s I + 2μE s , (6.15)

where  and μ are Lame constants. They are linked via shear modulus G, modulus of elasticity
(Young’s modulus) E and Poisson coefficient  as:
E E
μ = G= and  = . (6.16)
2(1 +  ) (1 +  )(1 2)
If the Hooke’s law is expressed in terms of the deviatoric stress (3.26) and strain tensors
(4.48):
1 1
 d =   (tr )I and E sd = E s  trEs I ,
3 3 ( ) (6.17)

78
 6. CONSTITUTIVE RELATIONS

one gets:
 2 
( ) 1
 d =   + μ trE s I + 2μE sd  (tr )I ,
3 3
(6.18)

and since tr d = trE d = 0, and trI = 3, it follows that:

 2 
( )
s 1
  + 3 μ trE  3 (tr ) = 0 . (6.19)

From (6.18) and (6.19) one gets:

 d = 2μE sd or  ijd = 2μijd
, (6.20)
( )
tr = ( 3 + 2μ ) trEs or p = 3Kv

( )
where p = 1 / 3(tr ) is pressure, v = 1 / 3 trE is the mean volumetric strain (see Section
s

4.27), and K is the bulk modulus given as:

2 E
K = + μ= , (6.21)
3 2(1 2)

6.3.2 Elastic potential and elastic constants’ limits

If thermal effects are neglected (divq=0, h=0), then, using (5.51) and (4.86), the energy
equation (5.46) can be written in the form:
u dEs
 = :D = : , (6.22)
t dt
and
s
DW = Du = :DE or DW = Du =  ij Dij . (6.23)

Assuming that W is a function of Es, we can write:

W
DW = Du =  ij Dij = D ij . (6.24)
ij

It follows from (6.24) that:

W
 ij = . (6.25)
ij

Function W which satisfies (6.25) is called elastic potential or strain energy function. Since the
reference (zero) level of strain energy function can be chosen arbitrarily, and the zero stress
state corresponds to the zero strain, the simplest form of strain energy function that leads to a
linear stress-strain relation is the quadratic form:
1
W = cijkl kl ij , (6.26)
2

79
 6. CONSTITUTIVE RELATIONS

which, by using (6.12), can be written as:

1 1
W = :Es or W =  ij  ij . (6.27)
2 2
The following principle results from the second law of thermodynamics: Among all states with
the same entropy, an isolated system will be in equilibrium in the state which possesses
minimum internal energy. Undeformed (or unstrained) state of ideal elastic solid is the state of
stable thermodynamic equilibrium (it will never spontaneously deform). The state of minimum
elastic potential has the corresponding (chosen) W=0, which means that W  0 must always be
satisfied.

Based on the expressions from the previous section, the strain-energy function (6.26) can be
expressed as:

( )
s 2 E  
( )

1 s s 2
W =  trE + μE :E =  trE s + Es :E . (6.28)
2 2(1+ ) 1 2 
As W must always be positive (vanishing only for zero strains), it follows that:
E

> 0 and > 0. (6.29)

2(1+ ) 1 2

or:
2
μ > 0 and  + μ = K > 0. (6.30)
3
(since we rule out the cases K=0 and μ=0). This implies that:
1
E > 0 and  1<
< . (6.31)
2
In practice, a material with negative  does not exist, while beryllium has value of  close to
zero. The value =0.5 implies that μ=E/3 and 1/K=0, or elastic incompressibility,
approximated by some rubbers.

6.3.3 Thermo-elasticity
In general, elastic coefficients in constitutive equation (6.12) are temperature dependent. For
small temperature variations this dependency is often neglected and the coefficients are
assumed constant. However, even then the thermal expansion of solids should be accounted for
if the resulting deformations are comparable to those due to mechanical loads, or if the thermal
deformations are constrained giving rise to thermal stresses. The thermal stresses can be
superimposed onto ‘mechanical’ stresses for linear materials (with linear constitutive
equations).

Constitutive relation which accounts for the temperature changes has the form:
 = C:E  B (T  T0 ) or  ij = cijkl kl   ij (T  T0 ) ,
s
(6.32)

80
 6. CONSTITUTIVE RELATIONS

where B is the second-order tensor, and T0 is the reference temperature corresponding to the
unstrained state of the material. Linear theory assumes that ij are independent of deformations
and that cijkl are independent of temperature. Then, for an isotropic solid at temperature T0 cijkl
must have the form (6.2), while for an isotropic material and  ij = 0 tensor B must have the
form I, and the constitutive relation (6.32) takes the so called Duhamel-Neumann form of the
Hooke’s law:

( s)
 =  trE I + 2μE   (T  T0 )I or
s
, (6.33)
 ij =  kkij + μ ij   (T  T0)ij

where  is the coefficient of volumetric thermal expansion, and is related to the coefficient of
linear thermal expansion  by the following relationship:
E
= = (2μ + 3 ) = 3K . (6.34)
1 2

6.4 Ideal and Newtonian fluids

6.4.1 Static and thermodynamic pressure

Experience indicates that a fluid at rest or in uniform flow cannot sustain shear stresses. Hence,
in these cases, shear stress magnitude is zero and the stress state is purely hydrostatic:
 =  p0I or  ij =  p0 ij , (6.35)

where p0 is the static pressure equal to the mean normal pressure defined by (3.25), with the
overbar added to distinguish from the thermodynamic pressure p, and with the subscript zero
to indicate that it corresponds to the fluid at rest or in uniform flow.
Thermodynamic pressure is defined by an equation of state called the kinetic equation of state
(to distinguish from the caloric equation of state):
p = p(,T ) , (6.36)

where  is the density and T is the absolute temperature. An example of an equation of state is
the ideal-gas law:
p = RT , (6.37)

where R is the gas constant. Notice that definition (6.36) does not guaranty equality between p
and p0 = 1 / 3(tr ) .

6.4.2 Ideal fluid

Ideal or inviscid fluid is defined as the fluid with shear stresses equal to zero even when it
flows. The constitutive relation for such fluid is:
 =  pI or  ij =  pij , (6.38)

81
 6. CONSTITUTIVE RELATIONS

where p is the thermodynamic pressure, which satisfies the state equation (6.36) [or more
precisely (6.37)]. No real fluid is inviscid, but the expression (6.38) has allowed a number of
analytical solutions for idealised flows.

6.4.3 Newtonian fluid

 +pI as a
The constitutive relation of a viscous fluid is expressed by so called viscous stress S=
function of the rate of deformation tensor D. If this relation is linear:
 =  pI + C:D or  ij =  pij + cijkl Dkl , (6.39)

such fluid is called Newtonian fluid. Since both the viscous stress tensor S and the rate of
deformation tenor D are symmetric, the constitutive relation for an isotropic fluid, similarly to
the isotropic solid, takes the form:
 =  pI +  (trD)I + 2μD or  ij =  pij + Dkk ij + 2μDij , (6.40)

where  and μ are two independent parameters which characterise fluid viscosity.
Constitutive equation of Newtonian fluids takes an interesting form when expressed in terms
of deviators:
1
 d =  + pI and Dd = D  (trD)I . (6.41)
3
Namely, expression (6.40) becomes:
 2 
 d = ( p  p)I +   + μ (trD)I + 2μDd , (6.42)
3

 d=0 and trDd=0, and trI=3, it follows from (6.42) that:

and since tr

( p  p) +  + 3 μ (trD) = 0.

2
(6.43)

Substituting (6.43) in (6.42), the first two terms on the right hand side vanish, and one obtains
the following two relations equivalent to (6.40):
 d = 2μD d
d
or  ij = 2μDij
1 D , (6.44)
p = p  trD = p  3 or p = p  Dkk = p + 3
 Dt

where p = 1 / 3(tr ) is the mean pressure, =1/3(trD) is the mean rate of volumetric dilatation
[see (4.93)], and
2
 = + μ, (6.45)
3
is the volumetric viscosity. The second form of the second equation in (6.44) employs:
1 D
trD = divv =  , (6.46)
 Dt

82
 6. CONSTITUTIVE RELATIONS

which follows from the continuity equation (5.17).

6.4.4 Stokes condition

Equation (6.43) shows that the mean pressure p is equal to the thermodynamic pressure p if
and only if one or both of the following conditions are satisfied:
1 D
trD = divv =  = 0 (Incompressible fluid),
 Dt
(6.47)
2
 =  + μ = 0 (Stokes condition).
3
The first condition is satisfied for incompressible fluids, where p is always equal to p. Strictly
speaking, there is no real fluids that are incompressible. However, in many practical situations,
the compressibility can be neglected for many fluid (flows). This is also the case for some
gasses whose speed of flow is much lower than the speed of sound. The second condition
 =  + 2 / 3μ = 0 which gives p = p , where  is the volume or bulk viscosity, is called Stokes
condition. It guaranties that pressure p can be defined as the mean normal stress even for
compressible fluids. This way, the thermodynamic pressure is expressed in terms of
mechanical stresses.
In fluid mechanics, almost without exception, it is assumed that the Stokes condition is
satisfied, i.e. =–2/3μ, and the constitutive equations (6.40) take the form:
2 2
 =  pI  μ(trD)I + 2μD or  ij =  pij  μDkkij + 2μDij
3 3 , (6.48)
d
 =  pI + 2μD or  ij =  pij + 2μDijd

which is equivalent to:

 d = 2μD d or  ij = 2μDijd
. (6.49)
p= p or p = p

The coefficient μ is called the dynamic viscosity. Generally, it is a function of temperature, and
in fluids it drops as the temperature rises, while in gasses it increases with temperature.

6.5 Fourier’s law and the caloric equation of state

The heat conduction in materials is defined by the relationship:

T
q =  gradT or qi =  ij , (6.50)
x j

where q is the heat flux, and  is a second-order tensor. In case of an isotropic material K=kI,
giving:

83
 6. CONSTITUTIVE RELATIONS

T
q = kgradT or qi = k , (6.51)
x i

where k is thermal conductivity. Expression (6.51) is known as Fourier’s law of heat

2
conduction, and is identical for solids and fluids. For the term 1 / T q  gradT to have non-
negative contribution to the entropy change [see (5.68)], k must be positive.
In cases where thermal effects cannot be neglected, caloric equation of state must be added to
the above presented (mechanical) constitutive relations:

( )
u = u T, ij , (6.52)

which gives the internal energy u as a function of T and ij for solids, and:
u = u(,T ) , (6.53)

giving the internal energy u as a function of  and T for fluids (see Section 6.26).

84
 7. MATHEMATICAL MODELS

Chapter 7

Mathematical models
7.1 Introduction
In this chapter, mathematical models, i.e. closed systems of equations, are formulated for
simple linear elastic solid and Newtonian fluid by combining the fundamental laws from
Chapter 5, which are valid for any continuum, and constitutive relations from Chapter 6 valid
for linear solids and fluids. Specifying initial (for unsteady problems) and boundary
conditions finally completes the models.

7.2 Linear elastic solids

First, basic equations are given for linear theory of elasticity employing balance equations
based on Lagrangian approach and Hooke’s law.

Basic equations

• Balance equations in integral form:

D
Dt 
dV = 0,
V
D
Dt 
vdV =   ndS +  fbdV , (7.1)
V S V
D
udV =   q  ndS +  hdV +  (:gradv)dV .
Dt 
V S V V

85
 7. MATHEMATICAL MODELS

• Balance equations in differential form:

D( J )
= 0,
Dt
Dv
 = div  + fb, (7.2)
Dt
Du
 = divq + h +  :gradv.
Dt

• Constitutive equations:
s
[( ) s
]
 = 2μE +  trE  ( 3 + 2μ ) (T  T0 ) I (thermo-elastic solid), (7.3)

q = kgradT (isotropic material). (7.4)

Hooke’s law is sometimes expressed in the following form, which is obtained by solving
(7.3) explicitly for Es:
1+  

Es =   (tr )  T  I . (7.5)
E  E 

• Equation of state:
u = cE T . (7.6)

Mathematical model
Hooke’s law assumes small strains so that the differences between material and spatial
descriptions can be neglected. Now, the basic equations of the linear theory of elasticity can
be written in the following forms:
• Integral form:

 t  t  dV =  {μ[(gradu) ] }
  u

+ gradu + [divu  ( 3 + 2μ) (T  T0 )]I  ndS

T

V S (7.7)
+  fbdV ,
V

 (c ET ) 
 dV =  kgradT  ndS +  hdV   (3 + 2μ)T0 (divu)dV . (7.8)
t t
V S V V

86
 7. MATHEMATICAL MODELS

• Differential form:
  u 
 {[ ]}
 = div μ (gradu) + gradu + grad[
divu  (3
+ 2μ ) (T  T0 )] + fb , (7.9)
t t
T

( c E T ) 
= div(kgradT ) + h  ( 3 + 2μ )T0 (divu) . (7.10)
t t

• Index notation:

  ui     u u j     u 
  = μ
i +  +  k  ( 3 + 2μ) (T  T0 )  + fb,i , (7.11)
t t x j   x j x i  xi x k 

( c E T )   T    u 
=  k  + h  ( 3 + 2μ)T0  k  . (7.12)
t x j  x j t  x k

It can be seen that the resulting system of 4 equations (in either integral or differential forms)
consists of 4 unknown functions of spatial coordinates and time ui(x,t) and T(x,t). Equation
(7.8) [or (7.10)] is called coupled energy equation, since its last term contains displacements
and therefore is coupled to the momentum equation. However, this term can be neglected in
most linear thermo-elastic problems. In such cases, the solution of a thermo-elastic problem
reduces to solving two separate (de-coupled) problems: i) solving of energy equation to
obtain the temperature field, and ii) solving of the momentum equations to obtain the
displacement field and, by using Hooke’s law, the stress field.
If , μ, , , k and cE have constant values, then the equations (7.9) and (7.10) can be written
in the form:
 2u
= μ u + ( + μ )grad(divu)  (3 + 2μ )gradT + fb ,
2
 (7.13)
t 2
T 
= k T + h  ( 3 + 2μ)T0 (divu) .
2
c E (7.14)
t t

87
 7. MATHEMATICAL MODELS

7.3 Newtonian fluids

Here, due to the nature of fluid flow basic equations are given for linear Newtonian fluids
employing balance equations based on Eulerian approach.

Basic equations
• Balance equations in integral form:

 t dV +  v ndS = 0,
V S
 (v )
 dV +  v( v n)dS =    ndS +  fbdV , (7.15)
t
V S S V
( u)
 dV +  uv ndS =   q  ndS +  hdV +  (:gradv)dV .
t
V S S V V

• Balance equations in differential form:

 D
+ div (v ) = 0 =  div v,
t Dt
( v) Dv
+ div( vv) = div  + fb or  = div  + fb, (7.16)
t Dt
(u ) Du
+ div (uv ) = divq + h +  :gradv  = divq + h +  :gradv.
t Dt

• Constitutive equations:
2
 =  pI  μdivvI + 2μD (Newtonian viscous fluid), (7.17)
3

q = kgradT (isotropic material).

• Equations of state:
Const. Incompressible fluid

 = (p,T ) = ( p) Barotropic fluid (7.18)

 p / RT Ideal gas

u = u(,T ).

88
 7. MATHEMATICAL MODELS

Mathematical model for Newtonian incompressible fluid

Using constitutive relation (7.17) and the equations of state for incompressible fluid
(=const., u=cvT), one gets the following mathematical model:

• Integral form:

 v ndS = 0 , (7.19)
S

 (v )

V
t
S S
{ [
dV +  vv  ndS =   pI + μ (gradv) + gradv ndS +  fbdV ,
T
]} V
(7.20)

( c v T )
 t
dV +  c v Tv  ndS =  kgradT  ndS +  hdV +
V S S V
(7.21)
 μ[(gradv) ]
T
+ gradv :gradvdV .
V

• Differential form:
divv = 0 , (7.22)

( v)
t {[
+ div( vv) = gradp + div μ (gradv) + gradv + fb ,
T
]} (7.23)

(c v T )
+ div (cv Tv) = div( kgradT ) + h +
t (7.24)
[
μ (gradv) + gradv :gradv .
T
]
• Index notation:
v j
=0, (7.25)
x j

v i (
 v iv j ) p   v
 μ i +
v j


t
+
x j
= +
x i x j   + fb,i , (7.26)

x j xi  


 (c v T )
+
( )
 cv Tv j
=
  T   v
+ h + μ i +
v j  v j
x j  k
x j  x j  . (7.27)
t  x j x i  x i

89
 7. MATHEMATICAL MODELS

One can see that the resulting system (in any of the forms presented above) consist of 5
equations with 5 unknown functions of spatial coordinates and time vi(x,t), p(x,t) and T(x,t).
It can also be noticed that for the incompressible fluids, the continuity (7.22) and momentum
(7.23) equations, which are also called the equations of hydrodynamics, are decoupled from
the energy equation (7.24).
Finally, if μ, k and cv have constant values, then the equations (7.22) to (7.24) can be written
in the following form:
divv = 0 Continuity equation (7.28)

Dv
 = gradp + μ 2 v + fb Navier-Stokes equation (7.29)
Dt

DT 2
c v = k T + h +  Heat equation (7.30)
Dt

where:

[ ] μ
[ ]
2
 = μ (gradv) + gradv :gradv = (gradv)T + gradv = 2μD:D,
T
(7.31)
2

is (always positive) dissipation function, and squaring of the rate of deformation tensor
should be understood as the double dot product.

7.4 Initial and Boundary Conditions

To complete the mathematical models, initial and boundary conditions have to be specified.
The number and type of initial and boundary conditions have to be selected to make sure that
the problem is mathematically well posed. A problem is considered to be well posed if its
solution exists and is unique, and depends continuously upon the initial and boundary
conditions. That is, a small perturbation of those conditions should give rise to a small
variation of the solution at any point of the solution domain and not be uncontrollably
amplified.

7.4.1 Classification of partial differential equations

The mathematical theory of partial differential equations (PDE) arising from the analysis of
problems of continuum (so-called equations of mathematical physics: wave equation, heat
conduction equation, Laplace equation,...) essentially deals with the way a particular

90
 7. MATHEMATICAL MODELS

perturbation propagates in time and space. By using a concept of characteristics (defined as

families of lines, surfaces or hypersurfaces, depending on the dimension of the problem,
along which certain properties remain constant or certain gradients become discontinuous),
we can distinguishes between three modes of propagation: hyperbolic, parabolic and elliptic.
In order to illustrate those concepts, let us consider a two-coordinate situation (time and one
space coordinate, or two space coordinates).

Hyperbolic problems

Hyperbolic PDEs arise from wave propagation phenomena. Information in a hyperbolic

problem travels at a finite speed called wave speed. An example is a rod with a force
suddenly applied at its end. An observer at a point along the rod will not be aware of it until
the wave reaches him. Discontinuities in initial and boundary conditions propagate through
the domain, and in non-linear hyperbolic PDEs discontinuities may develop even for
‘smooth’ initial and boundary conditions (e.g. shocks in compressible flow). A hyperbolic
problem has two real characteristics and the value of a dependent variable at point P (Fig.
7.1) depends on values of all points in region ABP (region of dependence), while the value in
P influences the values in the region PCD (zone of influence).

C zone of influence D
of point P

P charasteristics

2  3

region of
dependance
of point P
t A B
1
x

Figure 7.1 Region of dependence and zone of influence of point P for a hyperbolic problem

91
 7. MATHEMATICAL MODELS

Parabolic problem
Parabolic PDEs have solutions that are smooth in space (may possess singularities at
corners). An example of a parabolic equation is the heat conduction equation. Information
travels at an infinite speed in a parabolic system. For example, if heat source is applied at one
end of a rod, the temperature rises instantaneously along the entire rod according to the heat
conduction equation. Far from the source, the temperature increase may be very small. For a
parabolic problem two characteristics are identical and the region of dependence of point P
reduces to the line segment AB (Fig. 7.2). The zone of influence of point P, on the other hand,
extends to the whole region at one side of the characteristics AB.


2 3
zone of influence
of point P
A
B
P 1
t region of
dependance
of point P
x

Figure 7.2 Region of dependence and zone of influence of point P for a parabolic problem

Elliptic problem
Examples of elliptic PDEs are the Laplace equation. Here the solutions are smooth (even if
boundary conditions are ‘rough’). Also, boundary data at any point affect the entire solution.
The effect of small irregularities in boundary data tend to be confined to the boundary: this is
known as St Venant’s principle. Elliptic problems have no real characteristics and the
solution at P depends on all points in the solution domain , and inversely, the values of a
dependent variable at any other point in  is influenced by the value at P. In the other words,
the region of dependence is identical with the zone of influence and both of them coincide
with the whole solution domain  (Fig. 7.3).

92
 7. MATHEMATICAL MODELS



P

region of dependance 
y zone of influence 

Figure 7.3 Region of dependence and zone of influence of point P for an elliptic problem

As mentioned above, the classification of PDEs depends on whether lines or surfaces exist
across which the derivatives of the solution are discontinuous. This is equivalent to
examining whether lines exist along which PDEs can be reduced to ordinary differential
equations. Consider a second-order general linear partial differential equation in the domain x
and y:
2u 2u 2u
A( x,y ) + 2B( x,y ) + C ( x,y ) +
x 2 x y y 2
(7.32)
u u
a( x,y ) + b( x, y ) + c ( x, y )u = f (x,y).
x y

where u is the unknown variable. Depending on the equation discriminant:

2
D = B  AC , (7.33)

it can be shown that the following applies:

D > 0 – Hyperbolic problem with two families of characteristics and discontinuous derivatives,

D = 0 – Parabolic problem with one family of characteristics and smooth solution,

D < 0 – Elliptic problem with no real characteristic and with smooth solution.

7.4.2 Initial and boundary conditions

The steady momentum (for both Hookean solids and Newtonian fluids) and energy equations
are elliptic or, more appropriately, they are elliptic in space, requiring boundary conditions to

93
 7. MATHEMATICAL MODELS

be specified at all solution domain boundaries  . A range of boundary conditions is

applicable, but they can all be classified into two groups:
1. Dirichlet boundary conditions where the value of dependent variables (x), p(x) and T(x),
where (x) stands for u(x) and v(x), at the boundary is given (e.g. prescribed
displacement, velocity, pressure or temperature)
(x B ) = f1 , p(x B ) = f1 p, T (x B ) = f1T , x B D , (7.34)

2. Neumann boundary conditions where the value of a dependent variable gradient at the
boundary is specified (e.g. prescribed traction or given heat flux)
grad  (x B ) = f 2 , grad p(x B ) = f 2 p, grad T (x B ) = f 2T , x B N , (7.35)

where  D and  N are parts of the boundary  with prescribed either Dirichlet or
Neumann boundary conditions, respectively. In the case of the momentum equation for
solids, so-called mixed boundary conditions are applicable. Here, for example, the Dirichlet
boundary condition (given displacement) is applied in the direction normal to the boundary
surface, and the Neumann boundary condition (prescribed traction) is applied in the direction
tangential to the boundary surface.

In many practical situations the problem exhibits some kind of symmetry in which case only
a part of the body may be taken as the solution domain and the symmetry boundary
conditions are applied at the interface
 t ( x B ) p ( x B ) T ( x B )
n (x B ) = 0, = 0, = 0, = 0, x B  s , (7.36)
n n n
where n is direction normal to the symmetry plane  s , and n and t are displacement or
velocity vector components in the direction normal and tangential to the symmetry plane,
respectively.

The unsteady momentum equation for solids is hyperbolic or, more appropriately, it is
hyperbolic in time, and to start the calculations displacements and displacement rates at time
t = to have to be given at all points of the solution domain 
u (x,t0 )
u( x,t0 ) = u0 (x ), = u̇0 (x ), x  . (7.37)
t
The unsteady momentum equation for fluids and the unsteady energy equation are parabolic
or, more appropriately, they are parabolic in time, and require velocity and temperature to be
known at all points of the solution domain  at the initial instant of time t = to

v( x,t0 ) = v (x ), T ( x,t0 ) = T ( x), x  .

0 0
(7.38)

94
 7. MATHEMATICAL MODELS

7.5 Summary
It is important to note that equations governing momentum and energy transport in a
continuous media (see Sections 7.2 and 7.3 for thermo-elastic solids and Newtonian fluids,
respectively), can all be written in the form of a following generic transport equation:


t  
B dV +  B v ndS =   grad  ndS +  qS ndS +  qV dV , (7.39)
V S S S V

while the continuity equation is combined with momentum equation to obtain an equation for
pressure (in fluids usually). Here,  stands for the transported property, i.e. displacement u,
velocity v or temperature T. The meaning of the properties B and  are given in Table 7.1.
The term q S contains parts of the heat flux vector or the stress tensor, which are not
included in  grad , while q V contains the volumetric source terms.

Table 7.11: The meaning of various terms in the generic transport equation (7.39)
 B  q V q S
T cT k :gradv + h 0
u
μ(gradu) + [divu  ( 3 + 2μ ) (T  T0 )]I
T
u μ fb
t
 pI + μ (gradv)
T
v v μ fb

In case of thermo-elastic solid, the continuity equation does not have to considered (it serves
for determining density which is constant). In most practical applications it can be assumed
that not only the strains are small but also velocities, thus the (non-linear) convection term
vv , as well as the mechanical energy dissipation term  : gradv can be neglected.

95
 8. FINITE VOLUME DISCRETISATION

Chapter 8

Finite Volume Discretisation

In general, there is no method of obtaining the unknown functions
 =  (x,t ) and p = p(x,t ) (8.1)

satisfying eqs. (7.19) and (7.39) [or (7.19) and (7.7), (7.8), (7.20), (7.21)] in a closed form,
except in very few cases. Here,  stands for displacement u, velocity v and temperature T.
Therefore an approximate solution has to be sought (Fig. 8.1). If one assumes a certain variation
of dependent variables in time and space, the governing equations may be integrated and if it
happens that the assumed profiles are exact, these equations will be satisfied identically,
together with prescribed initial and boundary conditions. However, in general this is not the
case and the solution of the governing equations (7.19) and (7.39) amounts to guessing the exact
solution. It is not difficult to imagine that this is not an easy task. In order to make these profile
assumptions feasible the solution domain (time and space) is subdivided into a finite number of
time steps and control volumes, such that even a very simple profile assumption (e.g. linear)
will produce plausible results. Then a procedure (numerical method) is devised of utilising this
to obtain the unknown functions (8.1).

All numerical methods consist of transforming the governing equations into a system of (in
general non-linear) algebraic equations, with subsets of these approximating each conservation
equation. In order to achieve this, time, space and equations have to be discretised. The time
discretisation understands a subdivision of a given time interval into a number of smaller
subintervals or time steps. The space discretisation consists of defining a numerical mesh
consisting of a finite number of computational points, thus replacing the continuous functions
(8.1) by their values at those points. Replacement of individual terms in the governing equations
by algebraic expressions connecting nodal values on a numerical mesh is called equation
dicretisation.

It is intuitively expected that the accuracy of the numerical approximation depends directly on
the number of computational points, that is the better the discretised time and space represent

96
 8. FINITE VOLUME DISCRETISATION

the continuum, the better the numerical approximation. In other words, the error of a numerical
simulation has to tend to zero when the number of computational points tends to infinity
(discretised time and space approach the actual time and space), and if that is the case the
numerical discretisation is said to be consistent. The speed of this variation is characterised by
the order of the numerical discretisation.

MATH. MODEL ?? SOLUTION

(closed system =(x,t)
of equations) p=p(x,t)

CLOSED-FORM NUMERICAL
SOLUTION SOLUTION

Time Equations disc. Space

disc. (integ. & 'prof.' appr.) disc.
Dead end

Constant Variable Struct. Unstruct.

dt dt mesh mesh

Numerical considerations Physical implications

Accuracy Stability Efficiency Conservativeness Boundedness

consist. & implicit memory

order of & &
discret. explicit CPU time

Figure 8.1 Main concepts and ideas of numerical analysis

The structure of the algebraic equations resulting from the discretisation procedure leads to two
classes of methods, the explicit and the implicit methods. In explicit methods the unknown
variable can be explicitly expressed in terms of known values. However, this advantage is
seriously offset by the fact that stability of the method imposes severe restrictions on the
maximum admissible time step. While this might not be a limitation for physically unsteady
problems, where the small time step might be needed for the required time accuracy, it leads to
a large number of time steps in order to reach the steady-state solution corresponding to a

97
 8. FINITE VOLUME DISCRETISATION

physically steady problem. In implicit methods a solution of a system of equations is required,

since more than one set of variables are unknown at the same time level. In the majority of
cases, however, the structure of the coefficient matrix of that system is rather simple, such as a
block three-diagonal, five-diagonal etc. Thus, although the number of operations required will
be higher in implicit than in explicit methods, this is compensated by the fact that many implicit
methods (at least for linear problems) have no stability related time step restrictions, and hence
a smaller number of time steps will be needed to reach the steady state.

From the practical application point of view a very important property of a numerical method is
its efficiency i.e. computer memory and CPU time required to obtain a solution†. Although this
is influenced by the discretisation procedure, as well as by the physical properties of the system,
the choice of the solution algorithm for the system of algebraic equations resulting from the
discretisation is the most important step for the economy of a numerical method.
Apart from the above numerical aspects, one has to consider physical implications of the
discretisation scheme, like conservation of physically conserved quantities or boundedness of
the solution.

Some aspects of the Finite Volume (FV) method will be discussed in the next Section, where
1D energy and momentum equations are taken as model equations and the decisions regarding
the discretisation practice will be made. Based on those decisions, in Section 8.2 the full
derivation of the discretised equations for 2D situations will be conducted. Finally, FV
discretisation of 3D problems in domains of arbitrary shape will be outlined in Section 8.3.

8.1 Main Concepts and Ideas of the Finite Volume Method

8.1.1 Time and space discretisation

In order to obtain the discrete counterparts of the governing equations the time is discretised
into an arbitrary number of time steps of the size t and the body is subdivided into a finite
number of contiguous control volumes (CV) or cells of volume V bounded by the surface S
consisting of a number of cell faces of the area Sk (Fig. 8.2). The computational points (nodes)
are placed at the centre of each CV, while boundary nodes, needed for the specification of
boundary conditions, reside at the centre of boundary cell faces. No distinction is made between
the deformed and undeformed configurations since the strains are assumed to remain
sufficiently small.

† For example, Cramer's rule is a perfect mathematical model for the solution of the system of linear algebraic
equations. However, in order to solve a system of 30 equations, ~1033 arithmetic operations have to be performed,
which for the fastest present day computer exceeds by far the present estimate of the life of the Universe. Not to
mention that the solution would be useless due to an excessive round-off error propagation.

98
 8. FINITE VOLUME DISCRETISATION

computational node
boundary node

CV

Figure 8.2 Finite Volume space discretisation.

8.1.2 Equation discretisation

In order to illustrate the equation discretisation and problems associated with it, the 1D heat
conduction equation and 1D momentum equation will be considered, because of their simplicity
and easy physical interpretation.

Steady 1D heat conduction (problem definition)

We shall start with steady heat conduction in a rod (or in a channel filled with stationary fluid)
with fixed temperature at its ends (Fig. 8.3), i.e. we shall assume that the temperature is a
function of one space coordinate only
T = T ( x) , (8.2)

in which case the problem reduces to the solution of energy equation (7.8) [or (7.21)] reduced to
the following form
T
 k x n x dS +  hdV = 0. (8.3)
S V
Flux term Source term

where, in general, k = k(x), h = h(x) are known functions, with boundary conditions
T (0) = T0, T ( L) = TL . (8.4)

99
 8. FINITE VOLUME DISCRETISATION

T0 T(x) TL

x
L
Figure 8.3 Heat conduction in a rod (or channel)

Space discretisation

For simplicity, the solution domain is subdivided into N equal control volumes (Fig. 8.4) with
L
x i = = x (i = 1,2,...,N ) , (8.5)
N

xi
Sw Se

W P E
1 2 i-1 i i+1 N-1 N
xi

Figure 8.4 A rod (channel) subdivided into N equal control volumes

( Se = Sw = S, Vi = Sxi , n x,e = 1, n x,w = 1).

and the following notation is used

T ( xi ) = Ti = TP ,
T ( xi 1 ) = Ti 1 = TW , (8.6)
T ( xi +1) = Ti +1 = TE .

In that case eq. (8.3) becomes

T T
 k x dS   k x dS +  hdV = 0 . (8.7)
Se Sw V

Qe Qw HP

Calculation of integrals and gradients

The integrals in (8.7) are calculated by employing the Mid-Point rule and the following
approximation is obtained
 T   T 
ke   Se  kw   Sw + hPV = 0 . (8.8)
x e x w

100
 8. FINITE VOLUME DISCRETISATION

where the mid-point values of conductivity and heat source are, simply
 x   x 
ke = k xi +
  , kw = k xi   , hP = h (x i ) . (8.9)
2 2

To calculate temperature gradients in (8.8), the variation of temperature between the

computational points (the temperature profile) has to be assumed. If the simplest zero order
polynomial is assumed (Fig. 8.5a), it can be seen that the temperature gradient T / x at the CV
cell faces is not defined. If one supposes a piece-wise linear variation of temperature between
the grid points (Fig. 8.5b), then
 T  Ti +1  Ti TE  TP  T  Ti  Ti 1 TP  TW
 x = = ,  x = = . (8.10)
e x x w x x

a) b)

T
T

w e w e
i-1 i i+1 i-1 i i+1
x x
Figure 8.5 Temperature profile approximation: a) zero order polynomial, b) first order polynomial.

Discretised equations
By employing (8.10) eq. (8.8) becomes
ke k
x
( Ti +1  Ti )S  w (Ti  Ti 1 )S + hPVi = 0.
x
(8.11)

Qi = Qe Qi–1 = Qw Hi = HP

and if the following notation is introduced

S S
ai,i +1 = ke , ai,i 1 = k w ,
x x
 
S S
(
ai,i =   ke + kw  =  ai,i +1 + ai,i 1 ,
x x
) (8.12)

bi = H i =  hPVi

eq. (8.11) gets the form:

ai,i 1Ti 1 + ai ,iTi + ai,i +1Ti +1 = bi (i = 1,2,...,N ) . (8.13)

101
 8. FINITE VOLUME DISCRETISATION

By writing such an equation for all N control volumes a system of N algebraic equations with N
unknowns (T1, T2,..., TN) is obtained
a11T1 + a12T2 = b1
a21T1 + a22T2 + a23T3 = b2
a32T2 + a33T3 + a34 T4 = b3
......................................... .......
ai,i 1Ti 1 + ai ,iTi + ai,i +1Ti +1 (8.14)
= bi
......................................... .......
aN 1,N 2TN 2 + aN 1,N 1TN 1 + aN 1,N TN = bN 1
aN ,N 1TN 1 + aN ,N TN = bN

If conductivity k and heat source h (i.e. the coefficients) are not functions of T, eqs. (8.14) are
linear. If, however, either k or h is a function of T, eqs. (8.14) are non-linear.

The fact that for the calculation of temperature at any point three grid points are involved, which
renders this discretisation scheme implicit, is illustrated in Fig. 8.6, where the so-called
computational molecule is presented. Since maximum three non-zero coefficients are present in
anyone equation, system (8.14) is said to be sparsely populated or, simply, sparse, or to be more
precise three-diagonal.

1 2 i-1 i i+1 N

Figure 8.6 Computational molecule for steady state heat conduction equation.

Before we proceed to describe some methods for the solution of system of equations (8.13), and
to discuss some important issues linked to the steady heat conduction, namely conservativness
and boundedness, it is useful to write the equation of heat balance for a control volume (8.13) in
the following so-called compass notation form
aP TP = aETE + aW TW + bP , (8.15)
where
S S
aE = ke , aW = kw ,
x x
aP = aE + aW , (8.16)
bP = HP = hPVP .

102
 8. FINITE VOLUME DISCRETISATION

Solution of the system of discretised equations

If system of discretised equations (8.14) is linear, then it can be directly solved by employing a
number of linear equation solvers. Here, we shall mention only two methods. First is the
Thomas or three-diagonal matrix algorithm (TDMA), which is obtained by applying the
Gaussian elimination to the (three-diagonal) system (8.14), and is therefore a direct algorithm.
The second one is the Gauss-Siedel method (GSM), the best known and probably the most
efficient of point-iterative procedures for large systems of equations.
If system (8.14) is not linear, then it is usually solved by an iterative procedure, which includes
linearisation (assumption that coefficients and the right hand side vector components are
constant – calculated by using temperatures from the previous iteration), and then solution of
that system by one of the methods for systems of linear algebraic equations.

Let us, for clarity, write system (8.14) by employing the following notation
wi Ti 1 + di Ti + ei Ti +1 = bi (i = 1,2,...,N ) . (8.17)

TDMA uses the Gaussian elimination to put system (8.17) into an upper triangular form by
computing the new diagonal coefficients di and the new right hand side vector components bi as
w w
di = di  i ei 1, bi = bi  i bi 1 (i = 2,3,...,N ). (8.18)
di 1 di 1

Note that in the expression (8.18) the equals sign means "is replaced by". The values of
unknown temperatures are then calculated from the back substitution
b b eT
TN = N , Ti = i i i +1 (i = N  1,N  2,...,1). (8.19)
dN di

It can be seen that the TDMA requires both computer storage and time proportional to N, rather
than to N2 and N3, respectively, as it is generally the case with the Gaussian elimination. This
feature makes the TDMA a very powerful algorithm for solving large systems of linear
equations with a three-diagonal matrix.

The GSM consists of the following iterative procedure

Ti
(k+1)
=
1
(
di i
(k+1) (k)
b  wi Ti 1  ei Ti +1 ) (i = 1,2,...,N, k = 0,1,2,...) . (8.20)

where k is the iteration counter, and w1 and eN are taken as zero. The method is extremely
simple, but only converges under certain conditions. It is interesting to note that a sufficient
condition for the convergence is the diagonal dominance of the system coefficient matrix
( aP  aE + aW ) which will be shown shortly to be one of the requirements for the boundedness
of the solution.

103
 8. FINITE VOLUME DISCRETISATION

8.1.3 Conservativness
If the heat flux at each cell face is represented by the same expression for the two adjacent
control volumes, then the discretisation (numerical method) is conservative, because the fluxes
at the internal cell faces cancel each other and the overall heat balance is satisfied exactly
(except for round-off error) irrespective of the number of CVs (computational points) and not
only in the limit of the large number of computational points [see eq. (8.11)]:

– Qo+ Q1 + H1 = 0
– Q1+ Q2 + H2 = 0
– Q2+ Q3 + H3 = 0
...... .... (8.21)
–Qi–1+ Qi + Hi = 0
...... ....
– QN–1+ QN + HN = 0
---------------------------------------------------------
N
– Qo + QN +  Hi = 0
i =1
Although inherent in the FVM, the conservativness can be lost due to careless discretisation. If,
for example, in calculating the heat flux through the cell face between the cells P and E the cell-
centre instead the cell-face conductivity is used:
S
flux calculated from the cell P Qi = kP (Ti +1  Ti )
x
S
the same flux calculated from the cell E Qi 1 = kE (Ti +1  Ti )
x
will not cancel each other, since in general kP  kE, and the conservativness is lost.

8.1.4 Boundedness
Another important property of a discretisation scheme is its boundedness, which is closely
linked to the physical plausibility of the numerical solution. For example, in the absence of heat
sources (h = 0) no temperature should lie outside the range established by the boundary
temperatures, i.e. for To  TL:

T0  Ti  TL (i = 1,2,...,N ). (8.22)

Also, an increase in a boundary temperature (or temperature of any other point, with other
conditions remaining the same) must cause an increase (not decrease) of the adjacent point
temperature.

The discretisation formula (8.15) satisfies these requirements because all the coefficients have
the same sign (in this case positive). In addition to that, the central coefficient aP is equal to the

104
 8. FINITE VOLUME DISCRETISATION

sum of the neighbour coefficients aE and aW [eq. (8.16)] which means that if boundary
temperatures were increased by a constant, temperatures at all computational points would
increase by the same value. Also, equal boundary temperatures (To = TL) would produce a
uniform temperature field (Ti = To, i = 1,2,...,N).

There are, however, numerous formulations that may lead to coefficients of different signs. For,
example, if the source term HP is a function of temperature, then in order to enable the solution
of the system (8.14), it is often linearised, i.e. expressed in the following form
HP = C0 + C1TP . (8.23)
where Co and C1 are, in general, functions of temperature. In that case eqs. (8.15) and (8.16)
become
aP TP = aETE + aW TW + C0 ,
S S
aE = ke , aW = kw , (8.24)
x x
aP = aE + aW  C1 .

If C1 is negative, aP remains positive, i.e. of the same sign as aE and aW. However, if C1 is
positive, aP may become negative resulting in an unbounded solution (and in this case unstable
numerical algorithm). A positive C1 (in all points) clearly indicates a physically unstable
situation, because an increase in TP leads to the increase in the heat source, which in turn leads
to a larger TP, and so on.

Unbounded solutions are frequently encountered in cases of convection dominated flows. This
is illustrated in the next sub-section. Other situations of a potentially unstable and/or unbounded
numerical solution will be discussed in the Section on stability.

Steady 1D heat convection-conduction (problem definition)

Consider now a steady 1D flow of a conducting inviscid fluid through a channel of a cross
section S without heat sources (Fig. 8.7).

T(x), v(x)

. .
m m

x
L
Figure 8.7 Heat convection-conduction in channel

In this case the energy equation (7.21) reduces to the following form:

105
 8. FINITE VOLUME DISCRETISATION

T
 cTvn x dS =  k x n x dS . (8.25)
S S
Convective term Flux term

Assume that the fluid properties (

, c, k) are constant and that the velocity field v = v(x) is
known and it satisfies the continuity equation:

 vn x dS = 0 . (8.26)
S

Obviously, in case of a channel with constant cross sectional area, v=const. If the uniform space
discretisation (8.5) is used, then eqs. (8.25) and (8.26) become:
T T
 cTvdS   cTvdS =  k x dS   k x dS . (8.27)
Se Sw Se Sw

 vdS   vdS = (vS )e + ( vS ) w = ṁe + ṁw = 0. (8.28)

Se Sw

If the peace-wise linear variation of temperature between the grid points (Fig. 8.5b) is
employed, eq. (8.27) can be written as:
TP + TE T + TP k k
(cv )e 2
S  ( cv ) W
w 2
S = e (TE  TP )S  w (TP  TW )S
x x . (8.29)
Te Tw

By introducing the following notation

S c
aE = ke  ( vS )e ,
x 2
S c
aW = k w + (vS )w,
x 2
aP = ke
S
x
+ kw
S c
x 2 [
+ ( vS )e  ( vS ) w = ] (8.30)

[ ]
aE + aW + c ( vS )e  ( vS ) w = aE + aW
ṁe + ṁw = 0

the discretised equation becomes

aP TP = aETE + aW TW . (8.31)

Since the last term in the aP coefficient is equal to zero due to the continuity equation (8.28),
the property aP = aE + aW is preserved. However, the coefficients aE and aW are not
unconditionally positive, and whenever
S c
k   v S < 0,
x 2

106
 8. FINITE VOLUME DISCRETISATION

i.e.
v x
Pe = > 2, (8.32)
D

where Pe is the Peclet number and D = k / (c) is the thermal diffusivity, aE or aW become
negative and an unbounded solution may occur. This may also cause divergence of the Gauss-
Siedel solver (in the limit of zero conductivity the aP coefficients become zero!). For example,
3
consider air (
=1 kg/m , c=1000 J/kgK, k=0.025 W/mK) flowing with constant speed v=0.25
m/s through a channel with a cross section S=1 m2. If x=0.005 m, then Pe=50 and
aE = 120, aW = 130, aP = 10 .

and
a) for TE = 0 and TW = 0 TP = 1300
.
b) for TE = 100 and TW = 0 TP = 1200

3
which is unbounded and, obviously, unrealistic. Situation with flow of water (
=1000 kg/m ,
c=4200 J/kgK, k=0.5 W/mK) will be much worse, even for very small velocities.

a) The upwind scheme

A well known solution for the high Peclet number flows is to use the zero order polynomial
approximation (Fig. 8.5a) for the discretisation of the convection term in the so-called upwind
manner, e.g.
TP if flow is from P to E i.e. ṁe  0
Te =  , (8.33)
TE if flow is from E to P i.e. ṁe < 0

and similarly for Tw. To avoid the conditional statement in (8.33), the use of min and max
functions has proven useful, e.g.
ṁeTe = max(ṁe ,0)TP + min(ṁe ,0)TE , (8.34)

where the following relationship holds

max(ṁe ,0) + min(ṁe ,0) = ṁe . (8.35)

This results in the following form of discretised equation (8.27)

c[max(ṁe ,0)TP + min(ṁe ,0)TE ] +
k k (8.36)
c[max(ṁw ,0)TP + min(ṁw,0)TW ] = e (TE  TP )S  w (TP  TW )S,
x x

or by employing (8.35)

107
 8. FINITE VOLUME DISCRETISATION

 S S
ke + ke + c [ṁe  min(ṁe ,0)] + c[ ṁw  min(ṁw ,0)] TP =
  x  x

(8.37)
 S   S 
ke x  c min(ṁe ,0) TE +  k w x  c min(ṁw ,0) TW ,
   

which reduces to the form (8.31) if

S
aE = ke  c min(ṁe ,0),
x
S
aW = k w  c min(ṁw ,0), (8.38)
x
aP = aE + aW + c (ṁe + ṁw ) = aE + aW .

Thus, by using the lower order polynomial approximation of the cell-face temperature in the
convection term, a bounded discretisation ( aE > 0, aW > 0, aP = aE + aW ) is obtained.
However, this is achieved by sacrificing accuracy of the higher order approximation.

b) Blending of two schemes

To recover (some of) the lost accuracy and at the same time maintain the boundedness, a
number of schemes has been proposed, all of which amount to a blend of one lower order and
one higher order scheme, e.g.

Te = Te
LO
( )
+  e TeHO  TeLO , (8.39)

where 0  e  1 is the blending factor, which can either be fixed or the solution dependent. One
such scheme is the hybrid scheme, which takes:
TeHO if Pee  2
Te =  . (8.40)
TeLO if Pee > 2

Another scheme, which is preferred here, uses a constant blending factor to blend the first order
polynomial approximation with a small amount of the zero order one, e.g.
 T + TE 
ṁeTe = [max( ṁe ,0)TP + min( ṁe ,0)TE ] +  ṁe P  [max( ṁe ,0)TP + min( ṁe ,0)TE ] =
 2

[max(ṁe,0)TP + min(ṁe,0)TE ] +   ṁeTP + 2e (TE  TP )

ṁ
(8.41)

 ṁeTP + min(ṁe ,0)TP  min(ṁe ,0)TE ]

Implicit part (coefficient) Explicit part (right hand side vector)

and similarly for the west cell face. It is important to mention that the blending scheme is
implemented according to the ‘deferred-correction’ practice, i.e. convective flux is split into an
implicit part containing the contribution from the LO only (part of the coefficient), and an

108
 8. FINITE VOLUME DISCRETISATION

explicit part consisting of the difference between HO and LO contributions (part of the right
hand side vector), leading to the equation:

aP TP = aETE + aW TW + b , (8.42)

with

S
aE = ke  c min(ṁe ,0),
x
S
aW = k w  c min(ṁw ,0),
x (8.43)
aP = aE + aW + c (ṁe + ṁw ) = aE + aW ,
ṁ  ṁ 
b = c
min(ṁe ,0)  e (TE  TP ) + c
min(ṁw,0)  w (TW  TP )

2
2

One can see that the coefficients are unconditionally positive and that aP = aE + aW , which
enables the use of iterative solvers like Gauss-Siedel method, because the troublesome HO
contribution is relegated to the right hand side and can be controlled by choosing an appropriate
. However, this means that the solution procedure has to be iterative even if equation (8.42) is
linear.

8.1.5 Stability
A stable discretisation scheme is one for which errors from any source (round-off, truncation,
mistakes) are not permitted to grow in the sequence of a numerical procedure. Various
mathematical methods for stability analysis are available, the most popular being the Von
Neumann method based on a Fourier analysis. However, we shall limit ourselves to simple
physical considerations.

Unsteady 1D heat conduction (problem definition)

The concept of stability is in strict sense applicable only to marching (time dependent)
problems. For that reason the one-dimensional unsteady heat conduction will be analysed, i.e.
the temperature will now be considered to be a function of the spatial coordinate x and time t
T = T ( x,t ) , (8.44)

and the governing equation (7.8) [or (7.21)] will get the form
 T
t   k x n x dS  hdV .
cTdV = + (8.45)
V S V
Transient term (R) Flux term (Q) Source term (H)

109
 8. FINITE VOLUME DISCRETISATION

Here, in general, k = k(x,t) and h = h(x,t), with an initial condition

T ( x,t 0) = Tin ( x ) , (8.46)

and Dirichlet boundary conditions (see Fig. 8.2)

T (0,t ) = T0 (t ), T (L, t) = TL ( t) . (8.47)

Let us now discretise eq. (8.45) and examine the stability of the resulting equation.

Flux and source terms

We shall use the same spatial approximation for the terms on the right hand side of eq. (8.45) as
described in Section 8.1.2, keeping in mind that temperature is now a function of x and t:
T ( xi ,tm ) = Tim = TPm ,
T ( xi 1,tm ) = Ti 1 = TW , (i = 1,2,...,N; m = 1,2,...)
m m
(8.48)
T ( xi +1,tm ) = Tim m
+1 = TE .

i.e. we shall use the mid point rule to calculate integrals and assume a piece-wise linear spatial
variation of temperature (Fig. 8.4b) for the evaluation of gradients.

As for the temporal variation of the terms on the right side of eq. (8.45), the heat flux and the
source term will be assumed to be a weighted average of their values at the two successive time
levels, i.e. they will be evaluated at the time t = tm + (1  ) tm1 :
Qi = Qim + (1  )Qim1 = QPm + (1  )QPm1,
(i = 1,2,...,N; m = 1,2,...) (8.49)
Hi = H im + (1  )Him1 = H Pm + (1  ) HPm1,

where  is between 0 and 1.

Transient term
By employing the mid point rule the integral in the transient term is approximated as

 cTdV = (cT )i V = ( cT ) P V

m m
(i = 1,2,...,N; m = 1,2,...). (8.50)
V

For temporal variation we shall assume a two-time-level linear variation of temperature with
time and obtain the following approximation:

Ri =

t 
V
cTdV =
V
t
( [
cT )i  (cT )i
m m1
= ]
(8.51)

t [ ]
V
( cT ) m
P
 ( )
cT m1
P
(i = 1,2,...,N; m = 1,2,...).

110
 8. FINITE VOLUME DISCRETISATION

Discretised equations

Thus by using eqs. (8.49) and (8.51) and results from Section 8.1.2, the discretised form of the
eq. (8.45) can be written in the same form as eq. (8.13)
m m m m m m
ai,i 1Ti 1 + ai ,iTi + ai,i +1Ti +1 = bi (i = 1,2,...,N; m = 1,2,...). (8.52)

where

m = k m S S
ai,i +1 e , aim,i 1 = kwm ,
x x
 V m

ai,i =  (c )i + ai,i +1 + ai,i 1 ,

m m m
  t 
S S
ai,i +1 = (1  ) ke = (1  )k wm1 ,
m1 m1 m1
, ai,i 1 (8.53)
x x
m1 V
= ( c )i
m1
ai,i  aim1 m1
,i +1  ai,i 1,
t
[
m1 m1
bi =  ai,i m1 m1
1Ti 1 + ai ,i Ti ,i +1Ti +1 + H i + (1  )H i
+ aim1 m1 m m1
]
or in the alternative compass notation
m m m m m m m1 m1 m1 m1 m1 m1
aP TP = aE TE + aW TW + aE TE + aW TW + aP TP + bP . (8.54)
where
S m = k m S ,
aEm = kem , aW w
x x
V
m
aP =
t
( c ) mP + aEm + aWm ,
m1 S S
aE = (1  ) ke = (1  )k wm1 ,
m1 m1
, aW (8.55)
x x
V
aPm1 = ( c ) m1  am1 E  aW ,
m1
t P

bP = H Pm + (1  )H Pm1

Different values of temporal weighting factor  lead to different time discretisation schemes.
The best-known are the explicit scheme ( = 0), the Crank-Nicolson scheme ( = 0.5) and the
fully implicit scheme ( = 1), whose interpretation is given in Fig. 8.8. The only apparent
difference between the explicit and fully implicit schemes is that the former assumes that the old
m1 m1 m m
value QP and HP and the latter that the new value QP and HP prevail over the whole of
the time step. However, these two schemes have remarkably different properties. The Crank-
Nicolson scheme, on the other hand, assumes a linear variation of QP and HP with time and it
can be expected to be more accurate than the previous two.

111
 8. FINITE VOLUME DISCRETISATION

explicit
QPm-1= Qm-1
i

Crank-Nicolson
Q, H

Qm = Q m
P i

fully implicit

m-1 m
t

Figure 8.8 QP and HP variation for explicit, Crank-Nicolson and fully implicit schemes.

Explicit scheme
For  = 0 eqs. (8.54) and (8.55) can be written as
m m m1 m1 m1 m1 m1 m1 m1
aP TP = aE TE + aW TW + aP TP + H P . (8.56)
and
m1 m1 S m1 m1 S
aE = ke , aW = kw ,
x x
V
m1
aP =
t
( c ) m1 m1 m1
P  aE  aW , (8.57)
V
aPm =
t
( c ) m
P

One can immediately see that temperature at any computational point at the time tm is
dependent on temperatures from the previous time tm–1 only, which is illustrated by the
computational molecule in Fig. 8.9. This means that the coefficient matrix of the resulting
system of linear algebraic equations (obtained by linearisation of the non-linear equations, if the
original equations are non-linear) is diagonal, leading to a trivial matrix inversion and hence to a
solution with a minimum number of algebraic operations per time step and minimum computer
m
memory requirements (temperature TP can be directly calculated from eq. (8.56) and there is
no need to store coefficient matrix and right hand side vector). Moreover, if a steady state is
sought, one can still use the unsteady formulation (8.56) and march in time till the steady state is
reached. Any other scheme with   0 [as well as the steady state formulation (8.13)] would be
m m m
implicit, i.e. more than one temperature (e.g. Ti 1, Ti , Ti +1 ) is unknown at the same time and
the matrix to be inverted is not diagonal, hence requiring much more algebraic operations per
time step and much more computer memory.

112
 8. FINITE VOLUME DISCRETISATION

m+1
t

m-1
1 2 i-1 i i+1 N

Figure 8.9 Computational molecule for explicit scheme (o – unknown, • – known).

The question now arises: why do we bother with the other schemes? The answer is that we
wouldn't have if it was not for a serious limitation of the explicit scheme, which very often
m1
completely offsets its advantages. Namely, coefficient aP [eq. (8.57)] can become negative,
which can lead to an unbounded and unstable solution. For this coefficient to be positive, the
following condition (in case of uniform conductivity ke = kw = k, and V = S x) must be satisfied
at every instant of time
c(x) 2 k t 1
t < , or r = < . (8.58)
2k c (x) 2 2

This is the well-known stability criterion for the explicit scheme for parabolic differential
equations, which in many cases imposes severe restrictions on the maximum admissible time
step. Namely, the reduction of x required by the spatial accuracy causes maximum t to
decrease proportionally to x2, which makes the explicit scheme prohibitively expensive for
many practical applications. Since the stability of the explicit scheme depends on the parameter
r, it is classified as conditionally stable scheme.

To illustrate this, let us assume constant physical properties of a material and a situation without
sources and write eq. (8.56) in the following form
m
(
m1 m1
TP = r TE + TW )
+ (1 2r)TPm1 . (8.59)

Suppose that at the time tm–1: TEm1 = TWm1 = 100 K and TPm1 = 0 K (Fig. 8.10a). For r =1
formula (8.59) gives TPm = 200 K , which lies outside the bounds and is physically impossible.
m1
If temperature TP is increased to, for example, 200 K, one would expect to get an even
m m
higher temperature TP . However, eq. (8.59) produces TP = 0 K , which is again implausible. It
m1 m1
is interesting to note that for a uniform temperature field TE = TW = C eq. (8.59) gives the
m
correct answer TP = C .

113
 8. FINITE VOLUME DISCRETISATION

a) "Heater' b) 'Refrigerator'
200 0
t m t m
100 0 100 100 200 100
m-1 m-1

W P E W P E
x x
Figure 8.10 Physically unrealistic results produced by explicit time differencing scheme for r = 1.

Crank-Nicolson scheme
It can be shown that the Crank-Nicolson scheme (obtained for  = 0.5) is unconditionally
stable, but unbounded (physically unrealistic) solutions are possible. Namely, the coefficient
m1
aP (for uniform conductivity ke = kw = k, and V = Sx ) becomes [eq. (8.55)]
V  cx k 
m1
aP =
t
( c ) P  am1
m1
E
m1
 aW = S
t

x
(8.60)

m1
Thus, whenever, for a given space discretisation, the time step is not sufficiently small, aP
becomes negative and the solution may become unbounded. The stability of the scheme,
however, ensures that these unphysical oscillations eventually die out. It is important to note
that at any stage of the calculation 6 grid nodes are involved (3 at the current time tm and 3 at
the time tm–1), as illustrated by the computational molecule in Fig. 8.11. This in general requires
storing or recalculating of coefficients from the previous time level.

m+1
t
m

m-1
1 2 i-1 i i+1 N

Figure 8.11 Computational molecule for Crank-Nicolson scheme (o – unknown, • – known).

Fully implicit scheme
It can be observed from (8.55) that the only constant value of  which guarantees positive
coefficients is 1. Thus, for  = 1 eqs. (8.54) and (8.55) become

114
 8. FINITE VOLUME DISCRETISATION

aP TP = aETE + aW TW + bP , (8.61)
and
S S
aE = ke , aW = kw ,
x x
V
aP = aE + aW + (c) P , (8.62)
t
V
bP = (c )m1 TPm1 + HP
t P

where the superscript m denoting the new time interval has been omitted to facilitate the
comparison of eqs. (8.61) and (8.62) with the steady state equations (8.15) and (8.16).

Now it is easy to note one important practical feature of the fully implicit scheme. Namely, for
an infinite time step ( t   ) the fully implicit equations reduce to their steady state
counterparts. At the same time the fully implicit scheme involves only one previous time level
point, as can be seen from the computational molecule in Fig. 8.12, hence involving
significantly less computations and/or storage than the Crank-Nicolson scheme (not as evident
in this 1D case, as in 2D or 3D situations).

m+1
t

m-1
1 2 i-1 i i+1 N

Figure 8.12 Computational molecule for fully implicit scheme (o – unknown, • – known).
Additional considerations
Figure 8.13a shows the propagation of the solution outward from the initial conditions line in
the case of an explicit scheme. It can be observed that the value of the dependent variable
(temperature) at point P, three time steps away from the initial data line, is calculated without
any awareness of the boundary conditions. However, theoretical considerations in Section 7.4
indicate that the solution at point P should be influenced by the boundary conditions, since the
parabolic energy equation has the characteristics t = const.

115
 8. FINITE VOLUME DISCRETISATION

Thus, for a large t (or a small x, i.e. for r > 1/2) the solution is allowed to 'wander' (Fig. 8.13a)
and inevitably leads to unbounded physically unrealistic, unstable solution. For a smaller t (or
a larger x, i.e. for r < 1/2) the boundary conditions are able 'to catch-up with' the solution (Fig.
8.13b) and hold it inside the physical bounds, resulting in a stable procedure. On the other hand,
the implicit schemes involve boundary values in the calculation of temperature at all
computational points at every instant of time tm (m = 1,2,...) and therefore produce a stable
solution.

a)

P
m=3
t
m=2

m=1

m=0
i-1 i i+1
initial data
x line

b)

P
m=3
t
m=2

m=1

m=0
i-1 i i+1
initial data
x line

Figure 8.13 Propagation of the solution of explicit scheme.

116
 8. FINITE VOLUME DISCRETISATION

1D wave propagation (problem definition)

In the previous examples some properties of elliptic and parabolic problems were discussed. To
examine the properties of discretised hyperbolic equations, 1D stress wave propagation is
considered. Let us consider a bar of length L and constant cross section S made of a linear elastic
material with Young's modulus E. Horizontal movement of the left hand side of the bar is
constrained and a compressive stress pulse p of the duration tp is applied to the right hand side
of the bar (Fig. 8.14). The bar is assumed to be under uniaxial state of stress, i.e. only the axial
stress component xx exists and all other stress components are equal to zero.
u(x,t)
tp

p

x
L
Figure 8.14 1D bar loaded by compressive stress pulse

When a force is applied, its action is not transmitted instantly through all parts of the body. The
deformations produced by the force, and so the information about the instantaneous stress
distribution, travel through the body in form of stress waves. In the absence of thermal stresses,
Hooke's law [eq. (6.33) or (7.3)] takes the form
xx = E xx, (8.63)

and when body forces are omitted, the governing equation of motion (7.7) reduces to
 u u
t   dV =  E n x dS ,
t x (8.64)
V S
Inertial Surface
force ( Fin ) force ( Fx )

where u =u(x,t). The right hand side of eq. (8.64) is discretised in the same way as the first term
on the right hand side of heat conduction equation (8.45). To discretise the inertial part of eq.
(8.64), similar procedure to the one for the transient term of eq. (8.45) is used; i.e. by employing
the mid point rule the volume integral is approximated as
m
u  u   u
m
 u 
m
 t
 dV = 

 dV
t


= 
 t
i
V = 
 t V (i = 1,2,....,N; m = 1,2,....) (8.65)
P
V V i

and a two-time-level linear variation of the deformation velocity with time is adopted giving

117
 8. FINITE VOLUME DISCRETISATION

 u V   u m  u m1 
t  t
Fini =  dV =     =
t   t  i  t  i 
V
V
t
[( m
2  ui  ui
m1
) (
 ui
m1 m2
 ui = )]
(8.66)
V
t 2 [ uim  2uim1 + uim2 = ]
V m m1 m2
2 [uP  2uP + uP ] (i = 1,2,....,N; m = 1,2,....).
t
The surface force is evaluated at time t = tm + (1  ) tm1 as:
m m 1 m m  1.
Fxi = F xi + (1   )Fxi = Fx P + (1   )FxP (8.67)

Combining the equations (8.66) and (8.67), discretised counterpart of eq.(8.64) is obtained and is
of the same form as eq. (8.52) where temperature T is replaced by displacement u:
m m m m m m
ai,i 1ui 1 + ai,i ui + ai ,i +1ui +1 = bi (i = 1,2,....,N; m = 1,2,....) , (8.68)

where
m m S m m S
ai,i +1 = Ee , ai,i 1 = E w ,
x x
m
V 
ai,i =    2 + aim,i +1 + aim,i 1 ,
 t
m1 m1 S m1 m S
ai,i +1 = (1  )E e , ai,i 1 = (1  )E w ,
x x
V (8.69)
m1
ai,i = 2 2  aim1 m1
,i +1  ai,i 1,
t
 m1 m1 + a m1 um1   V um2

bi =   aim1 m1
,i 1ui 1 + ai ,i ui i,i +1 i +1 i .
 t 2 

Equation (8.68) can be written in the compass notation as:

m m m m m m m1 m1 m1 m1 m1 m1
aP uP = aE uE + aW uW + aE uE + aW uW + aP uP + bP , (8.70)
where
S S
aEm = E em , m
aW = E wm ,
x x
V
aPm =  2 + aEm + aW m
,
t
S m S
aEm1 = (1  )E em1 , m1
aW = (1  )E w ,
x x (8.71)
V
aPm1 = 2 2  am1 m1
E  aW ,
t
V m2
bP =  2 uP .
t

118
 8. FINITE VOLUME DISCRETISATION

In case of the explicit discretisation scheme (=0), the stability criterion ( aPm1 > 0) leads to a
following condition on the time step size:
x E t
t < , or r = < 1. (8.72)
E /  x

Here, E /  defines the speed of the wave propagation along the bar made of a linear elastic
isotropic material (Ee = Ew = E).

The computational molecule for fully implicit scheme is shown in Fig. 8.15. It involves only one
point from each previous time level (m–1, m–2) point, and therefore requires much less
computations and/or storage than the Crank-Nicolson scheme.

m+1
t

m-1

m-2
1 2 i-1 i i+1 N
x

Figure 8.15 Computational molecule for fully implicit scheme (o – unknown, • – known).
8.1.6 Accuracy
In Sections 8.1.2, 8.14 and 8.1.5 both surface and volume integrals are calculated by employing
the Mid-Point rule. In Section 8.1.2 the spatial temperature variation was approximated by a
piecewise-linear function of the spatial coordinate x leading to the expressions (8.10) for the
temperature gradient. Similarly, in Section 8.1.5 the time derivative was approximated by a two-
time-level approximation, which for different values of  led to the explicit, the Crank-Nicolson
and the fully implicit time differencing schemes. Finally, in Section 8.14 the cell-face values of
the temperature were obtained by blending first and zero order polynomial approximations.

Intergrals
Let us first establish the accuracy of the Mid-Point rule. If function f(x) is continuous and has at
least three continuous derivatives in [xw, xe], then it can be shown that

119
 8. FINITE VOLUME DISCRETISATION
xe
xe  x w n x 2 ''
 f (x)dx =
n  mi
f + ( xe  x w ) 24 f ( ),  ( xe ,x w ) , (8.73)
xw i =1
Exact value Mid-Point formula Error

where n is the number of subintervals, and fmi is the value of f(x) in the middle of the ith
subinterval. Since the error of the above formula is proportional to x2 [x = (xe – xw)/n], the
Mid-Point rule is second order accurate.

Cell-face gradients
In order to examine the accuracy of gradient approximations consider a continuous, sufficiently
times differentiable function f(x) and its Taylor series expansion (TSE) around the point xi

 f  1  2 f  2 1  f 
3
f (x) = f (xi ) +   ( x  x i ) +  2 ( x  xi ) +  3 ( x  x i ) + ...,
3
x i 2  x  6  x  (8.74)
i i
Then the values of f(x) at points xi+1 = xi + x and xi–1 = xi – x can be written as

 f  x + 1   f  x 2 + 1   f  x 3 + ...
2 3
f (x i +1 ) = f (xi ) +
 x
i 2  x 2
6  x 3

i i
(8.75)
 f  1  2 f  1 3 f 
f (x i 1) = f (xi )  x + 2
x 2  3
x 3 + ...
 x i 2  x 6  x
i i

By using equations (8.75) different approximations to (f/x)i can be obtained. For example,
from the first equation (8.75) it follows that

 f  f (xi +1)  f (x i ) 1  2 f  f  f
 x =  2
x + ... = i +1 i + O(x ) (8.76)
i x 2  x x
i
Exact Taylor series Truncation TSA TE
value approximation error

while the second equation (8.75) gives

 f  f (xi )  f (xi 1 ) 1  2 f  f  f i 1
 x = + 2
x + ... = i + O(x ) (8.77)
i x 2  x x
i
Exact Taylor series Truncation TSA TE
value approximation error

or by subtracting the second from the first equation (8.75) one gets
f (xi +1)  f (x i 1) 1  3 f  2
 f 
 x
i
=
2x 6  x
f  f
 3
x + ... = i +1 i 1 + O x
2x
2
( ) (8.78)
i
Exact Taylor series Truncation TSA TE
value approximation error

120
 8. FINITE VOLUME DISCRETISATION

Equations (8.76) and (8.77) correspond to the so-called forward and backward differencing
formulae, respectively, and are known under the common name one-sided approximations (or
one-sided differencing schemes). They are first order accurate, because their truncation error
(TE) is proportional to x. The last approximation [eq. (8.78)] corresponds to the so-called
central differencing scheme, and is second order accurate because its TE is proportional to x2.
It is interesting to note that, although formula (8.78) uses only two points, as the first-order
formulae (8.76) and (8.77), one gets second order accuracy. This is due to the fact that the
tangent at the point xi to a second order polynomial (parabola) fitted through points xi–1, xi, xi+1
is parallel to the secant through the end points xi–1 and xi+1 (Fig. 8.16). If, however, xi  (xi–1 +
xi+1)/2, the second order accuracy is lost.

tangent

f
secant
parabola

f i-1 (f i) f i (fe) fi+1

i-1 (i) i (e) i+1

2x (x)
x
Figure 8.16 Geometrical interpretation of the second order accuracy of the two-point formula.

This leads to the idea that formulae (8.76) or (8.77) can be regarded as the central differencing
formulae with respect to point xe = xi+1/2 = (xi + xi+1)/2 (values in brackets in Fig. 8.16) or xw =
xi–1/2 = (xi –1+ xi)/2 respectively, giving second order accuracy at these points.

Accuracy of the present discretisation

Thus, according to the preceding analysis, our spatial discretisation is second order accurate (on
a uniform mesh), while the temporal accuracy depends on the time discretisation scheme
--------------------------------------------------------
Scheme Truncation error
--------------------------------------------------------
Explicit scheme O(x2,t)
Crank-Nicolson scheme O(x2,t2)
Fully implicit scheme O(x2,t)
--------------------------------------------------------

121
 8. FINITE VOLUME DISCRETISATION

In fact, integrals and dependent variable gradients could be obtained by involving any number
of adjacent points in space and time with the order of approximation increasing with the number
of points. However, there is always a trade-off between the order of accuracy and the number of
points simultaneously (implicitly) involved in the calculation, which determines the bandwidth
of the system of linear algebraic equations to be solved to obtain Ti or ui (i = 1,2,...,N). In the
present equation discretisation, in case of a steady state or implicit and Crank-Nicolson
schemes, in the calculation of dependent variable at any computational point three grid points
are involved (Figs. 8.6, 8.11, 8.12, 8.15) resulting in a three-diagonal system of linear algebraic
equations of the form (8.13).

Dissipation and dispersion (Truncation error)

Although the truncation error tends to zero as the number of computational points and time
steps tends to infinity, for the finite number of computational points and time steps the
truncation error is always present and manifests itself in the form of what is known as
dissipation and dispersion. While dissipation (sometimes called diffusion) tends to smooth-out
sharp gradients, dispersion usually expresses itself in the form of non-physical oscillations (Fig.
8.17). In general, if the leading term in the truncation error contains an even derivative [as in
eqs. (8.76) and (8.77)], the solution will predominantly exhibit dissipative errors (Fig. 8.17b).
On the other hand, if the lowest order term in the truncation error is an odd derivative [as in eq.
(8.78)], the solution will predominantly exhibit dispersive errors (Fig. 8.17c).
a) b) c)

Figure 8.17 Effects of dissipation and dispersion on the computation of sharp gradients.
(a) Exact solution.
st
(b) Numerical solution distorted primarily by dissipation errors (typical of the 1 order methods).
nd
(c) Numerical solution distorted primarily by dispersion errors (typical of the 2 order methods).

8.1.7 Consistency
Consistency deals with the extent to which the discretised equations approximate governing
conservation equations. A discrete counterpart of a conservation law is said to be consistent if
the difference between the two vanishes as the time and space discretisation become sufficiently
fine (i.e. if t  0, x  0).

122
 8. FINITE VOLUME DISCRETISATION

Since the approximation involved in the present discretisation procedure is in approximating

integrals (Mid-Point rule), and the time and space gradients, and since the truncation error of
those approximations vanishes under the time step and mesh refinement, the present
discretisation is consistent.

8.1.8 Convergence
In the context of numerical continuum mechanics the term convergence is given various
meanings. In theoretical considerations convergence is usually related to the numerical method
on the whole which is said to be convergent if the solution of the discretised equations
approaches the exact solution of the original governing equations when t and x tend to zero,
i.e. when the number of computational points tends to infinity. It can be proved that for a well-
posed (linear initial value) problem and a consistent discretisation scheme, stability is the
necessary and sufficient condition for convergence.
In the present course the term convergence will be related to two iterative processes:
• (outer) iterations set up to solve the problem of non-linearity and/or inter-equation coupling,
and
• solver (inner) iterations performed during iterative solution of the system of linear algebraic
equations resulting from disretisation and subsequent linearisation and/or decoupling.

8.1.9 Summary
Let us now summarise the findings from the previous sections and make decisions about the
choices regarding the present FV method.
1. In order to achieve a conservative discretisation the diffusion fluxes have to be uniquely
associated with cell faces, and not cells.
2. The sufficient condition for a bounded (physically plausible) solution (in the absence of
sources) is that all the coefficients of the discretised equations are positive and that the
diagonal coefficient is equal to the sum of all the neighbouring (in both space and time)
coefficients.
3. Integrals are calculated by employing the Mid-Point rule, and the linear spatial variation of
dependent variables is assumed guaranteeing a second order accurate spatial discretisation.
4. The fully implicit temporal variation scheme is preferred on the grounds of its unconditional
stability, boundedness and better efficiency than, for example, the Crank-Nicolson scheme.
5. The truncation error of the adopted discretisation O(x2,t) vanishes proportionally to x2
(on a uniform mesh) and to t, which means that the present discretisation is consistent.

123
 8. FINITE VOLUME DISCRETISATION

8.2 Two-dimensional FV Discretisation

8.2.1 Two-dimensional Approximations

Thermo-elastic solid

In some practical situations the geometry and the (mechanical and thermal) loading
configuration allows reasonably accurate results to be obtained by reducing the problem to a
two-dimensional one. These situations are known as (Fig. 8.18):

a) plane strain problem,

b) plane stress problem,
c) axi-symmetric problem.
z
 r

z y
x
y
x
z

a) Plain strain b) Plain stress c) Axial symmetry

Figure 8.18 Two-dimensional approximations for solids.

In all these cases it is sufficient to consider a two-dimensional domain in x–y plane. In case of
axi-symmetric problems, radial direction is represented by x axis and axial one by y axis, and
corresponding displacement components by u and v, respectively. The momentum equations
[second equation in (7.1)] take the form
 u  
t  t
( )
 dV =   xx n x +  xy n y dS +  fbx dV  <  x
dV >,
V S V V
(8.79)
 v
t  ( )
 dV =   yx n x +  yy n y dS +  fby dV ,
t
V S V
with the stress tensor components
u  u v u

 xx = 2μ +  + + < >   [ 3KT ],

x  x y x 
v  u v u

 yy = 2μ +  + + < >   [ 3KT ], (8.80)

y  x y x 

124
 8. FINITE VOLUME DISCRETISATION

u  u v u

<   = 2μ +  + +  [ 3KT ] >,

x
x y x 
 u v

 xy =  yx = μ + ,

y x 
where the terms enclosed within < > brackets are included only for axi-symmetric cases, and
those within [ ] brackets account for thermal effects.

The energy equation [third equation in (7.1)] can be expressed as


t 
( )
cTdV =  q x n x + qy ny dS +  hdV , (8.81)
V S V

and the heat flux components

T T
qx = k , qy = k . (8.82)
x y

Now, these three situations will be discussed in some detail.

Plane strain problem

Situations involving a long prismatic body whose geometry and loading conditions do not vary
significantly in the longitudinal z direction is referred to as a (generalised) plane strain problem
(Fig. 8.18a). In these cases the dependent variables (displacement components and temperature)
can be assumed to be functions of only x and y coordinates and time t
u = u(x, y,t), v = v(x,y,t), T = T (x,y,t). (8.83)

If a further assumption is made that the displacement w in z direction is zero at every cross
section
w = w(x, y,z, t) = 0, (8.84)
then the strain components xz , yz and zz vanish and equations (6.33) and (7.7) reduce to
equations (8.79) and (8.80).
From equations (6.33) follows that xz = yz = 0, but the normal stress in the z direction is not
zero and can be easily calculated from (6.33) from the condition zz = 0

( )
 zz =   xx +  yy  [ ET ] (8.85)

Plane stress problem

In contrast to the plane strain condition, where the longitudinal dimension in z direction is large
compared to the dimensions in x and y directions, cases involving very small thickness and no

125
 8. FINITE VOLUME DISCRETISATION

loading in z direction are known as a (generalised) plane stress problem (Fig. 8.18b). It is,
normally, assumed that
 zz =  xz =  yz = 0 (8.86)

and the non-zero stress components xx, yy and xy are averaged over the thickness and
assumed to be independent of z. It can be shown that, after some algebraic manipulation, eqs.
(8.80) remain valid if and 3K are replaced by
E E E E
=  ; 3K =  . (8.87)
(1 +  )(1 2 ) 1  2 1 2 1 

It can be noted that the above formulae can formally be obtained by replacing 1 – 2 with 1 – 
in respective expressions for  and 3K in eqs. (6.16) and (6.21). However, although the strain
components xz and yz vanish, the component zz is given by expressions
 1+  
zz =  ( +  yy +
1  xx
)
1 
( )
T =   xx +  yy + T ,
E
(8.88)

obtained from (6.33) and the condition zz = 0.

Axi-symmetric problem
Many engineering problems are concerned with solids of revolution which are symmetrical in
shape about axial axis z and are subjected to axially symmetrical loading (Fig. 8.18c). It is
convenient to express these problems in terms of cylindrical coordinates r,
and z, shown in
Fig. 8.18c. Due to axial symmetry
u = 0,  r =  z = 0 and  r =  z = 0, (8.89)

where u
is circumferential displacement component. The non-zero strain components are
ur u u 1  u u 
rr = ,  zz = z ,  = r ,  rz = zr = r + z
,, (8.90)
r z r 
2 z r 

where ur and uz are displacements in r and z directions respectively. The non-zero stress
components are rr, zz,

and rz. Using eqs. (6.33) and (8.90) they can be expressed as
ur  u u u

 rr = 2μ +   r + z + r  [ 3KT ],
r
r z r
uz  u u u

 zz = 2μ +   r + z + r  [ 3KT ],
z
r z r
(8.91)
u  u u u

 = 2μ r +   r + z + r  [ 3KT ],
r
r z r
 u u

 rz =  zr = μ r + z ,

z r 

126
 8. FINITE VOLUME DISCRETISATION

Note that if r is replaced by x, z by y, and ur , uz by u, v, respectively, the above equations

reduce to equations (8.80).

Mathematical model for two-dimensional thermo-elastic stress analysis

By combining (8.79) to (8.82) the following system of equations is obtained (for 2D
approximations, the mathematical model of a thermo-elastic linear solid [eqs. (7.7) and (7.8)],
reduces to):
• Momentum
 u  u
u v u  
u v  
  dV =   2μ +   + + < >  [3KT ] n x + μ  +  n y dS +
t t  x x y x  y x  
V S  
1 u
u v u 
 fbx dV  <  x  2μ x +   x + y + x   [3KT ]dV >,
V V
(8.92)
 v 
u v   v
u v u   
  dV =  μ +  n x + 2μ +   + + < >   [ 3KT ] n y dS +
t t  y x  y x y x   
V S

 fby dV ,
V
• Energy

T T 
t  cTdV =   k
x
nx + k n y  dS +  hdV .
y 
V S V

Eqs. (8.92) form a closed system of 3 equations with 3 unknown functions u, v, T of spatial
coordinates x, y and time t, which makes a mathematical model for 2D thermo-elastic stress
analysis. Note that if the heat source term h in energy equation does not contain the thermo-
elastic dissipation term, the energy equation is decoupled from momentum equations.

Newtonian fluid
There are also situations when fluid flow can be considered to be two-dimensional. These are
a) planar 2D flow (e.g. flow in a rectangular duct, Fig. 8.19 a), and
b) axi-symmetric 2D flow (e.g. flow in a pipe, Fig. 8. 19b).

127
 8. FINITE VOLUME DISCRETISATION

a) b)

. . . .
m m1 m1+m 2

.
m2
Fig
ure 8.19 Two-dimensional flow in a) rectangular duct, and b) pipe.

In case of a 2D planar and axi-symmetric fluid flows, the mathematical model for an
incompressible Newtonian fluid [eqs. (7.19) to (7.21)] takes the following form:
• Continuity

 ( un x + vn y )dS = 0,
S
• Momentum
  u
u v 

t  ( )
udV +  u un x + vn y dS =   2μ n x + μ  +  n y  dS 
V S  x S
y x  

 pn x dS +  fb x dV ,
S V
(8.93)
 
u v  v

t  ( )
vdV +  v un x + uny dS =  μ  +  n x + 2μ n y dS 
y x  y 
V S S
v
 pn y dS +  fb y dV + <  2μ y 2 dV >,
S V V
• Energy

T T 
t  ( )
cTdV +  cT un x + vn y dS =   k
x
nx + k n y  dS +
y 
V S S
  

u 2
v 

u v 
2 2 2

v
μ 2  + 
   x   y    y x 
+ + + < 2   > dV +  hdV .
y
V     V

This is a closed system of 4 equations with 4 unknown functions u, v, p and T of spatial

coordinates x and y and time t. Note that terms in <> brackets correspond to axi-symmetric
situations and are omitted in case of a planar flow.

128
 8. FINITE VOLUME DISCRETISATION

8.2.2 Discretisation for a 2D Cartesian Mesh

For reasons of simplicity and clarity the 2D equations from the previous Section are discretised
in this Section. The discretisation of the full three-dimensional equations would only bring in
much longer expressions whose sheer size would obscure the underlying methodology.

Time and space discretisation

The observed time interval is divided into a number of equal sub-intervals or time steps t and a
two-dimensional solution domain which can be sub-divided into N control volumes by a non-
uniform Cartesian numerical mesh (e.g. Fig. 8.20) is considered. In that case the (outward)
surface unit vectors have the coordinates (Fig. 8.21)
1   1  0 0 
ne =
, nw =
, nn =
, ns =
 . (8.94)
 0 0  1   1

Figure 8.20 Examples of 2D solution domain

nn
Sn

n
Se
nw
w CV e
ne
y Sw
s

Ss ns
x

Figure 8.21 Control volume and control surface S =  Sk (k = e,w,n,s) .

Equation discretisation

According to the findings from Section 8.1, for the temporal variation a fully implicit temporal
differencing scheme is used, i.e. all dependent variable values are related to the current time tm,
except for ones coming from the time derivative term, which are denoted by an appropriate
superscript (m–1 for the time tm–t, and m–2 for the time tm–2t). For the spatial variation a

129
 8. FINITE VOLUME DISCRETISATION

central differencing scheme is used to express the cell-face values of gradients of dependent
variables in terms of their neighbouring nodal values. For the cell-face values of the dependent
variables featuring in the convection terms the SO scheme blended with a small amount of the
FO scheme will be implemented in a deferred-correction manner. All integrals are calculated by
employing the Mid-Point rule.
Different from the practice employed in Section 8.1, a non-uniform grid spacing will be used
here (Fig. 8.22), implying that the spatial interpolation factor of 0.5 used for calculation of the
cell-face values of dependent variables and physical properties from the two neighbouring nodal
values will be replaced by the following spatial interpolation factors:
Pe Pw Pn Ps
fe = , fw = , fn = , fs = . (8.95)
PE PW PN PS

xn = xs = xP

NW NE
N
nw n ne
yn

y e = y w =yP
W w P e E

sw s se y s

S
SW SE
x w x e

Figure 8.22 Two-dimensional control volume and compass labelling scheme.

Thermo-elastic solid
a) Momentum
a.1) u-momentum:
Now, the equilibrium of forces acting in the x direction on the CV shown in Fig. 8.23a [first
equation (8.79)] can be written in the following form
 u  
t  t
( )
 dV =    xx n x +  xy n y dS +  fb x dV  <  x
dV >, (k = e,w,n,s). (8.96)
V k Sk V V

Inertial Surface Body <Axi-symmetric

force force force contribution>

130
 8. FINITE VOLUME DISCRETISATION

By employing (8.80), eq. (8.96) becomes

 u  u  v u
  u v

t   dV =
t  (2μ + ) x +  y + < x >  [ 3KT ]dS +  μ y + x  dS 
V Se Sn
Fxin Fxe Fxw (8.97)
 u  v u
  u v

 (2μ + ) x +   y + < x
>   [ 3KT ]dS 
  μ  y + x  dS +
Sw   Ss
Fxn Fxs
1 u  u v

 fbx dV  <  x (2μ + ) x +  x + y   [3KT ]dV > .
V V
Fbx < Fx >

a) b)
Fyn Qn
Fxn
n
Fby Fye
Fxw Qw
w e
F bx Fxe Qe
Fyw CV CV
s

Fxs
F ys Qs

Figure 8.23 Forces a) and heat fluxes b) acting on a control volume.

Note that eq. (8.97) is exact, i.e. no approximation has been introduced so far! Before applying
the discretisation practices from Section 8.1 to approximate the individual terms in the above
equation, it should be noticed that apart from the approximation of the so-called normal
derivatives u/x at the east and west cell faces and u/y at the north and south cell faces, one
has to approximate the so-called cross derivatives v/y at the east and west cell faces and v/x
at the north and south cell faces, for which no counterpart exists in the energy equation. This is
done, for the north cell face for example, by interpolating between the gradients at points P and
N (Fig. 8.22) in the following way
 v   v   v   v  vw   v  vw 
 x  = (1 f )  x  + fn   = (1 f n ) e + fn  e , (8.98)
n
n
P x N x  P x  N

where
ve = (1 f e )v P + fe vE ,
v w = (1 f w )v P + f wvW .

131
 8. FINITE VOLUME DISCRETISATION

By employing the above and the discretisation practices described in Section 8.1, the forces
acting on the CV in the x direction can be approximated by:

Fxin 
V
t 2
(
uP  2uPm1 + uPm2 , )
S  v  ue
Fxe  (2μ +  ) e (uE  uP ) + Se
 + <  Se >  [3K Te Se ],
x e y  e xe
S  v  u
Fxw  (2μ +  ) w (uP  uW ) + Sw
 + <  w Sw > [ 3K TwSw ],
x w y  w xw
Sn  v 
Fxn  μ (u  uP )
y n N
+ μS n  ,
x n
(8.99)

S  v 
Fxs  μ s (uP  uS ) + μS s ,
y s x  s
Fbx  fbxP V,
1  uP  1 1  
< Fx  (2μ +  ) +  ( ue  u w ) + ( v n  v s )  [ 3K TP ]V > .
x P  xP  x P y P  
where the superscript m is the time step counter and the compass labelling scheme shown in Fig.
8.22 is used. The material properties μe, e etc. are obtained by linear interpolation between the
two neighbouring computational points, e.g.:
μe = (1 fe )μP + fe μE .

Finally, after introducing the following notation:

S S
aE = (2μ +  ) e ; aN = μ n ;
xe y n
S S
aW = (2μ +  ) w ; aS = μ s ;
xw y s
v  v  u u
b = Se    Sw   + <  e Se   w Sw > +
 y  e
 y  xe
w
xw
v  v 
μSn    μSs   (8.100)
x n x s

[ 3K Te Se  3K TwSw ] + f bxPV +  V2 (2um1

P  uP ) 
m2 
 t 
V  1 1  
<   ( ue  uw ) + ( v n  v s )  [3K TP ] >;
xP 
x P y P  
 V  V
aP =  aK +  2 + < (2μ +  ) 2 > (K = E,W ,N,S),
K  t  xP

the momentum equation (8.97) becomes

132
 8. FINITE VOLUME DISCRETISATION

aP uP   aK uK = b (K = E,W ,N,S ), (8.101)

K

Note that if thermal stresses are not present (T = 0), then expression in [] brackets is zero, if a
steady state is considered the expressions in {} brackets vanish and therefore are not calculated,
while terms in brackets <> are non-zero only in axi-symmetric cases.

a.2) v–momentum:
Similarly, from the equilibrium of forces acting in the y direction one gets:
S S
aE = μ e ; aN = (2μ +  ) n ;
xe y n
S S
aW = μ w ; aS = (2μ +  ) s ;
x w y s
 u
 u

b = μSe    μSw   +
y  e y  w
(8.102)
 u
 u
u us
Sn  S s + <  n Sn   S >
x  n x  s xn xs s

[ 3K Tn Sn  3K TsSs ] + f byPV +  V2 (2v Pm1  v Pm2 );

 t 
 V 
aP =  aK +  2  (K = E,W ,N,S ),
K  t 

and
aP v P   aK v K = b (K = E,W ,N,S ), (8.103)
K

b) Energy
Discretisation of the energy equation is even simpler. For the CV in Fig. 8.23b eq. (8.81)
becomes

t 
cTdV =   (q x n x + q y n y )dS +  hdV (k = e,w,n,s), (8.104)
V k Sk V

Rate of change Surface Heat

of internal energy heat flux source

or by using (8.94) and (8.82) (see Fig. 8.23b)

 T T T T
t   k x dS   k x dS +  k x dS   k x dS +  hdV
cTdV = (k = e,w,n,s). (8.105)
V Se Sw Sn Ss V
R Qe Qw Qn Qs H

133
 8. FINITE VOLUME DISCRETISATION

By employing the equation discretisation described in Section 8.1, one obtains

R
V
t[( cT )P  (cT )m1
P ]
,

S S
Qe  ke e (TE  TP ), Qw  k w w (TP  TW ),
x e x w (8.106)
S S
Qn  k n n (TN  TP ), Qs  k s s (TP  TS ),
y n y s
H  hPV,
and after introducing the following notation
S S
aE = ke e ; aN = k n n ;
x e y n
S S
aW = k w w ; aS = ks s ;
x w y s
(8.107)
V m1 
b = hPV +  ( cT )P ;
  t
V 
aP =  aK +  ( c ) P  (K = E,W ,N,S ),
K  t

the energy equation (8.105) gets the same form as momentum eqs. (8.101) and (8.103):
aP TP   aK TK = b (K = E,W ,N,S ), (8.108)
K

c) Initial and Boundary Conditions

To start the calculation, all dependent variables have to be set to their initial values: uo, vo and
To and in the case of the momentum equations the displacement vector components u1 and v1
at time t1 = tot are also required.
The expressions for the evaluation of forces and heat fluxes described above are valid for all
interior cell faces, which are shared by the two adjacent cells. Boundary cell faces, however, are
related to only one cell and require special treatment, which depends on the type of the
boundary conditions.

c1) Dirichlet boundary conditions

If, for example, an east cell face of a CV happens to be at the solution domain boundary and the
displacement is known (Fig. 8.24a), the uE and vE in the expressions (8.99) and in the
corresponding expressions for y direction are simply replaced by uB and vB. Similarly, if the
temperature TB at the east boundary is given (Fig. 8.24b), then TE in the expression (8.106) for
the heat flux Qe is replaced by TB.

Note that xe now stands for the distance between the points P and B, xB.

134
 8. FINITE VOLUME DISCRETISATION

a) b)
tBy , vB
t B, uB

P P qB
B=e tBx , uB B=e
TB

xB xB

Figure 8.24 Dirichlet and Neumann boundary conditions for momentum a) and energy b).

c2) Neumann boundary conditions

If, again, an east cell face is a boundary cell face and the traction tBx is supplied (Fig. 8.24a),
then expression for Fxe in eq. (8.99), and the analogous expression for y momentum Fye have to
be replaced by
Fxe = tBx Se , Fye = tBy Se . (8.109)

The values of the displacements uE and vE required for the calculation of gradients (v / y)e ,
(u / y)e etc. are obtained by extrapolation from the interior of the solution domain.
Similarly, if the heat flux qB is given (Fig. 8.24b), the Qe given by eq. (8.106) is to be replaced
by
Qe = qB Se . (8.110)

c3) Mixed boundary conditions

If, for example, at the east-face boundary (Fig. 8.25a) displacement uB in the x direction and
traction tBy in the y direction are given, then the u-equation is treated as having the Dirichlet
boundary conditions, and the v-equation is treated as having the Neumann boundary conditions,
i.e. uE in the expression for Fxe in (8.99) is replaced by uB, while the expression for Fye in the v-
momentum equation is replaced by tBySe. The value of displacement vE required for the
calculation of gradient (v / y)e is obtained by extrapolation.

135
 8. FINITE VOLUME DISCRETISATION

a) b) vB
t Bx
t By
B=n yB
P
B=e uB
P

xB

Figure 8.25 Mixed boundary conditions for momentum equations.

c4) Symmetry plane

If the problem considered is symmetric in such a way that the south boundary, for example, can
be taken for a symmetry plane, then by looking at a near-boundary cell and its mirror image on
the other side of the symmetry plane (Fig. 8.26), it can be seen that the following holds:

n
vP
w P e
TP uP
yB
Symmetry sw s=B se
plane
TS
S
uS
vS

Figure 8.26 Symmetry plane boundary conditions.

uS = uP , v S = vP ,
us = uB = uP , vs = v B = 0,
 u   u  v 
= , = 0,
 x  s  x  P  x  s (8.111)
and
TS = TP, y s = 2yB ,
Ts = TB = TP,

and the south boundary surface force components are:

136
 8. FINITE VOLUME DISCRETISATION

S  v

Fxs  μ s (uP  uS ) + μSs = 0,

y s
x  s
S  u
u
Fys  (2μ +  ) s (v P  v S ) + Ss
 + <  s Ss > [ 3K TsSs ] = (8.112)
y s x s xs
Ss S u
( 2μ +  ) v P +  s (ue  uw ) + <  P Ss > [ 3K TP Ss],
y B x s xs

while the heat flux is equal to zero

Ss
Qs  ks (T  TS ) = 0.
y s P
(8.113)

Newtonian fluid

Firstly, the convection term will be discretised for the case of 2D planar flow and a Cartesian
control volume shown in Fig. 8.21. After that discretised equations for energy and momentum
will be assembled.

a) Convection term discretisation

The convection flux of the variable  through the CV face k represents the rate at which variable
 is convected into (out of) the control volume through the cell face Sk (Fig. 8.21). Taking the
east cell face as an example, it may be approximated as follows:

  (un x + vn y ) dS  ṁee
*
Ce = , (8.114)
Se

where e* stands for the cell-face mean value of velocity components u, v or internal energy cT,
and ṁe is the mass flux through the east cell face

 u dS  ue Se .
*
ṁe = (8.115)
Se

The way cell-face values of velocity components ue* and ve* and convected variable e* are
calculated has a strong influence on both the accuracy and stability of a numerical method. The
* *
velocity components ue and ve are obtained from a special interpolation practice which assures
*
a stable solution procedure as described in the next Section, while e is calculated by blending
the second-order accurate, symmetric central differencing formula

 1 
eCD = (1 f e ) P + fe E = ( P + E ), for a uniform mesh fe = 0.5 (8.116)
 2 

with some small amount of the first order asymmetric (upwind) differencing scheme, which
assumes that the value of  at the cell face is calculated as

137
 8. FINITE VOLUME DISCRETISATION

P if flow is from P to E

UD
e = (see Section 8.14).; (8.117)
E if flow is from E to P

Hence

(
e* = eUD +   eCD  UD
e , ) (8.118)

where  is a blending factor with a value between zero and one, thus combining the accuracy of
a second order scheme and the stability of a first order one. It is important to mention that the
CD scheme is implemented according to the 'deferred-correction' practice, i.e. convective flux is
split into an implicit part containing the contribution from the UD only (part of the coefficient)
and an explicit part consisting from the difference between CD and UD contributions (part of the
right hand side vector, see Section 8.14).

b) Discretised equations
After assembling all the terms featuring in eqs. (8.93), they can be written in the form of the
following non-linear algebraic equation which is identical to the one derived for solid body
stress analysis (4.1)
aP P   aKK = b (K = E,W ,N,S), (8.119)
K

where  stands for u, v or T.

•b1) momentum

The coefficients for both u and v velocity components are the same
S S
aE = e e  min(ṁe,0); aN =  n n  min(ṁn ,0);
x e yn
S S
aW =  w w  min(ṁw,0); aS = s s  min( ṁs ,0); (8.120)
x w ys
 V

aP =  aK +   (K = E,W ,N,S ),
K t 

and the source terms are

bu =  u [min(ṁe ,0)  ṁe f e ] ( uE  uP ) +  u [min(ṁw,0)  ṁw f w ](uW  uP ) +
 u [min(ṁn ,0)  ṁn f n ](uN  uP ) +  u [min(ṁs,0)  ṁs f s](uS  uP ) +
S S S S
e e ( uE  uP )   w w ( uP  uW ) + n n ( uN  uP )  s s (uP  uS ) +
x e x w y n y s
V m1

n n (v ne  vnw )  s s (vse  v sw )  ( pe Se  pwSw ) + fbxP V + 

S S
uP 
x n x s t 
(8.121)

138
 8. FINITE VOLUME DISCRETISATION

bv =  v [min(ṁe ,0)  ṁe f e ] (v E  v P ) +  v [min(ṁw ,0)  ṁw f w ](vW  v P ) +

 v [min(ṁn ,0)  ṁn f n ](v N  v P ) +  v [min(ṁs,0)  ṁs f s](v S  v P ) +
S S S S
n n (v N  v P )  s s (vP  v S ) + e e (v E  v P )  w w (v P  vW ) +
y n y s xe x w
V m1

e e ( une  use )   w w (unw  usw )  ( pn Sn  psSs ) + fbyP V + 

S S
vP 
y e y w t 

b2) energy

The coefficients for the energy equation differ from the ones for the heat conduction in solids
(8.107) by the presence of the convection terms only
S S
aE = ke e  c e min(ṁe ,0); aN = k n n  c n min(ṁn ,0);
x e yn
S S
aW = k w w  c w min(ṁw,0); aS = ks s  c s min(ṁs ,0);
x w ys
V 
aP =  aK +  cP  (K = E,W ,N,S),
K  t
(8.122)
b = T ce [min(ṁe ,0)  ṁe fe ](TE  TP ) +  T cw [min(ṁw,0)  ṁw f w](TW  TP ) +
T c n [min(ṁn ,0)  ṁn f n ](TN  TP ) +  T c s[min(ṁs ,0)  ṁs f s ](TS  TP ) +
 V (cT ) m1

hPV + (:gradv) +  P 
t 

b3) continuity
In assembling equations (8.119) the continuity equation [see eq.(8.93)] in the discretised form
 ṁk = 0 (k = e,n,w,s), (8.123)
k
where the mass fluxes ṁk defined by (8.115) are multiplied by the value of dependent variable
(u, v, cT) at point P, and subtracted from the corresponding transport equation. In this way the
central coefficient remains to be the sum of the neighbouring coefficients (in both space and
time). An independent use of the continuity equation, i.e. its conversion in an equation for
pressure will be dealt with in the next Section.

c) Calculation of pressure. SIMPLE algorithm

It can be noted that in the above procedure pressure, featuring in the source terms of the
momentum equations (8.121), has remained unknown, while at the same time no independent
use has been made of the continuity equation (8.123). The problem lies in the fact that pressure
does not feature explicitly in the continuity equation which consequently cannot be considered

139
 8. FINITE VOLUME DISCRETISATION

as 'an equation for pressure' and it just comes as an additional constraint on the velocity field.
This constraint can be satisfied only by adjusting the pressure field. However, pressure is not a
conserved property and has no its governing transport equation, so it is not immediately clear
how this adjustment of pressure is to be performed.

At the same time, the pressure source term in the momentum equation is calculated in a way
that amounts to using a second-order space-centred scheme for the calculation of pressure
gradient [e.g. (8.78)]. For example, the pressure force acting on a control volume in the x
direction is calculated as (see Fig. 8.23a)
Fpx = (  pdS   pdS) = ( pe Se  pw Sw ) =
Se Se
(8.124)
1 1  1
  ( pP + pE )Se  ( pW + pP )Sw =  ( pE  pW )S (if Se = Sw = S)
2 2  2
where a uniform mesh spacing is assumed (fe = 0.5). One can see that only alternate (and not
adjacent) computational nodes are involved in the pressure force calculation and an oscillatory
(checkerboard) pressure field is interpreted as uniform by the momentum equation (Fig. 8.27).
Such a scheme can produce a correct pressure-difference field, although the underlying pressure
field possesses an unphysical oscillatory profile.

p2

1/2(p1+p 2)
p
p1

i-2 i-1 w i e i+1 i+2

x
Figure 8.27 A wavy pressure variation 'seen' as uniform by a second order scheme.

c1) Cell-face velocity

The simple and yet efficient way of getting round both these problems is to calculate the fluid
velocity at a cell face in the following manner, e.g. for the east cell face

*  V   p  pw pE  pW 
ue = ue +  v  e   , (8.125)
 aP  e  x P xe
v
where the overbar denotes a spatially interpolated quantity [formula (8.116)], aP the
corresponding momentum equation central coefficient and xP = xn = xs (Fig. 8.22). The
second term in eq. (8.125) is a third-order pressure diffusion term, whose role is to smooth out

140
 8. FINITE VOLUME DISCRETISATION

the oscillatory pressure velocity profile, and at the same time introduce pressure into the
continuity equation in such a manner that a pressure-correction equation can easily be
constructed. This is achieved by employing the predictor-corrector procedure defined by the
SIMPLE algorithm which will now be briefly described.

c2) Predictor stage; pressure-correction equation

The so called predictor stage values of u, v and p, featuring in expressions for the cell-face
velocity components like (8.125), which satisfy the (linearised) momentum equations, do not
necessarily satisfy the continuity equation (8.123). By correcting the cell-face velocities, e.g.
ue** = ue* + ue' (8.126)

by expressing the velocity corrections in terms of the pressure corrections by using a truncated
form of eq. (8.125)

'  V  p  pW
' '
ue =  v E , (8.127)
 aP  e xe

and requiring that the double-starred cell-face velocities satisfy the continuity equation (8.123),
an equation of the form (8.119) for pressure correction p' is obtained, with coefficients
V  S  V  Sn
aE =   v e , aN =   v ,
 aP  e x e  aP  n x n
(8.128)
 V  Sw V S
aW =   v , aS =  v s ,
 aP  w x w  aP  s  x s

aP =  aK (K = E, N, W , S) ,
K
b =   ṁk (k = e, n, w, s) ,
k
where all variables have their predictor stage values, i.e. those that satisfy the (linearised)
momentum equations.

c3) Corrector stage

After the field of pressure correction p' is obtained, velocity components and pressure are
corrected via
1 '
('
uP = uP,pred + v pe  pw Se ,
aP
)
aP
1 '
('
v P = v P,pred + v pn  ps Sn , ) (8.129)

p = ppred + pp',

141
 8. FINITE VOLUME DISCRETISATION

where p is an under-relaxation factor (typically p = 0.2 to 0.3), and is necessary because the
approximations introduced in deriving the pressure-correction equation often result in
overestimating the magnitude of p', which in turn leads to slow convergence or divergence of
the solution procedure. Finally, the mass balance satisfying fluxes are calculated as follows, e.g.
'
( '
ṁe = ṁe,pred + ṁe = ṁe,pred  aE pE  pP
'
) (8.130)

and used in the consequent calculations of convective fluxes.

c4) Closure
It should be noted that the right hand side of the pressure-correction equation is the mass
imbalance [see eqs. (8.128)], caused by the fact that the predictor stage velocity components do
not in general satisfy the continuity equation at the initial stage of calculation. When eventually
the continuity equation is satisfied, b is equal to zero at all computational points, which results
in a zero p' field and a converged solution is obtained. In fact, the p' field is of no concern
whatsoever, and these quantities only serve in a temporary way to procure a satisfactory linkage
between the pressure and velocity fields so that they ultimately satisfy both the momentum and
continuity equations. Thus, taking the truncated form of the expression (8.125) (or any other
approximation) during the derivation of the pressure correction equation is justified as long as
the procedure converges.

The nature and significance of the third-order pressure diffusion term in expression (8.125) are
identical to those introduced by expression (8.146) in Section 8.4.1 where the diffusion
transport through an arbitrarily oriented cell face is discussed. Once a converged solution is
obtained, this correction term will vanish if the pressure varies linearly or quadratically in space,
or it will become small in the limit of a fine grid. In that case the velocity at the cell face will be
a simple spatial average, implying that the cell-centre velocities will also satisfy the continuity
equation.

It is important to note that the coefficient matrix of the pressure correction equation is
symmetric, but singular. However, in case of a well posed problem (the global mass balance
satisfied exactly) the resulting set of linear algebraic equations is undetermined, which
practically means that the pressure correction, and consequently the pressure level is
undetermined. This, however, has no consequences on the solution, since only pressure
gradients (differences), and not the pressure itself, feature in fluid momentum equations [see
e.g. (8.124)].

d) Initial and boundary conditions

To start the calculation in an unsteady case all dependent variables have to be prescribed the
appropriate initial values.

142
 8. FINITE VOLUME DISCRETISATION

A number of different types of boundary conditions is encountered in fluid flow problems. The
most frequently used ones will now be discussed.

d1) Inlet
This is an inflow boundary where the distribution of all variables is know, and it is therefore of
a Dirichlet type.

d2) Outlet
This is an outflow boundary where the distribution of dependent variables is not known and a
kind of Neumann boundary condition is applied. Namely, the velocity components are
extrapolated in a way that forces an exact (to the round-off error) overall mass balance, e.g. if
outlet is at the east boundary
uB = uP +u, (8.131)
where

ṁinf   uP dS
SB
u = , (8.132)
S B

where ṁinf is the total mass inflow and SB is the total outlet area. The other dependent variables
are obtained by simple (usually zero order) extrapolation.

d3) Wall
For the momentum equation Dirichlet boundary conditions are usually applied, i.e. the boundary
velocity components are set equal to the wall velocity (zero in the case of a stationary wall),
while the boundary conditions for the energy equation can be either of Dirichlet (specified
temperature) or Neumann (specified heat flux) type.

d4) Symmetry plane

At planes of symmetry the 'mirror image' approach described previously is employed.

d5) Pressure-correction boundary conditions

Since the pressure correction is not a physical but numerical quantity, the pressure-correction
equation boundary conditions are derived from the physical boundary conditions. Thus, if
velocity at a boundary is known (as is the case in all the above boundary types), the velocity
correction is zero and from the definition of the velocity correction (8.127) it follows that the
pressure correction gradient is zero, implying Neumann type boundary conditions.

e) Turbulent flow
Observations of fluid flow have revealed that two different flow regimes, laminar and turbulent,
can be distinguished. While in a laminar flow the migration of fluid elements from one layer to

143
 8. FINITE VOLUME DISCRETISATION

another does not take place, i.e. the flow can be thought of as a motion in essentially parallel
paths, turbulent flow is characterised by a chaotic motion of fluid elements and their migration
not only to the neighbouring, but also to distant regions and it is inherently three-dimensional
and unsteady. It is also observed that a stable laminar flow occurs only at a very low range of
Reynolds numbers (in pipes Re =
uD/ < 2300) and otherwise is turbulent.

Stokes' constitutive relationship is valid irrespective of the flow regime. However, the small
length and time scales of turbulent flow make its direct solution absolutely impractical, and
some kind of turbulence modelling is often used to overcome that difficulty. The most
frequently employed is the so-called two equation k– model of turbulence, which solves a
transport equation of the form [3rd expression in (7.15) or in (8.93)] for both the turbulent
kinetic energy k and its dissipation rate , and then calculates the so-called turbulent viscosity
tur and turbulent conductivity ktur via simple algebraic relationships involving k and . The
governing transport equations (7.15 or (8.93)) retain their form, but u, v, p and T now stand for
the respective time averaged quantities and  and k for their effective values (eff =  + tur,
keff = k + ktur).

In conjunction with the k– model the so-called wall functions are used in formulating boundary
conditions at wall boundaries.

144
 8. FINITE VOLUME DISCRETISATION

8.3 Solution procedure

The present discretisation practice generally leads to a system of non-linear algebraic equations,
with sub-sets for each dependent variable, which has to be solved at each time step (Fig. 8.28).
If equations are non-linear, an iterative procedure has to be used. This frequently includes some
kind of linearisation of the equations, thus reducing the original problem to the series of linear
equation problems.

The inter-equation coupling (e.g. v displacement appearing in the u-momentum equation, and
vice-versa) is often another reason for an iterative solution procedure. Namely, instead of
solving a system of equations containing all dependent variables, the decoupled sub-systems for
each dependent variable are obtained by temporarily assuming all other dependent variables to
be known (taken from an initial guess or previous iteration time step). This approach, known as
segregated solution algorithm, is employed by the present FV method and will be presented in
the next Section. A part of that algorithm is the solution of the systems of linear algebraic
equations and the Incomplete Cholesky Conjugate Gradient (ICCG) is used.

8.3.1 An Iterative Segregated Procedure

Since the Finite Volume methods favour segregated over simultaneous solutions of discretised
algebraic equations, in this Section a segregated algorithm is described.
As a result of FV discretisation of conservation equations, and by integrating over a time
interval t and over all N CVs three in general mutually coupled sets of N non-linear algebraic
equations with three unknowns (u, v, T) of the form (8.119) are obtained for either thermo-
elastic solid or Newtonian fluid

aP P   aKK = b (K = E,W ,N,S),

K

where  stands for u, v or T. An iterative procedure for the solution of these equations will now
be presented.

145
 8. FINITE VOLUME DISCRETISATION

System of 3 x N
(non)linear alg.
equations

'Genuine'
iterative
methods

With decoupling (Linearisation) (Linearisation)

Without decoupl.
3 syst. of N eqn. With decoupling Without decoupl.
1 syst. of 3xN eqn.
3 syst. of N eqn. 1 syst. of 3xN eqn.

outer iterations
no
CONVERGED iterations Iterative Direct
inner

yes

no yes no
SOLUTION Converged CONVERGED

yes

SOLUTION

Figure 8.28 Various options for solving systems of discretised equations

• First, all dependent variables are given the corresponding initial values (actual initial values
in case of a transient analysis, and an arbitrary initial guess if only the steady state is sought).
Then the boundary conditions which correspond to the first time step are applied and the sets of
equations for each individual dependent variable are linearised and temporarily decoupled by
assuming that coefficients and source terms are known (calculated by using dependent variable
values from the previous iteration or the previous time step) resulting in a system of linear
algebraic equations of the form
A = b (8.133)

for each dependent variable, where A is an N x N matrix, vector  contains values of dependent
variable  at N nodal points and b is the source vector.

• These systems are then solved sequentially in turn until a converged solution is obtained.
The procedure is assumed converged when the following conditions are satisfied for all three
sets of equations

146
 8. FINITE VOLUME DISCRETISATION

N N
  aK K + b  aP P < p  aP P ,
i =1 K i =1 (8.134)
i( n)   (i n 1) < q i(n 1) (i = 1,2,...,N)

where p and q are typically of the order of 10–3 and superscripts (n) and (n–1) denote values at
the two successive iterations.

• In the next time step the whole procedure is repeated, except that the initial values are
replaced by the values from the previous time level.

Systems of linear equations [eq. (8.133)] can be solved by a number of methods. Convergence
of those so called inner (solver) iterations is promoted by diagonal dominance of the coefficient
matrix. Since (coefficients and) sources are only approximated (based on the values of
dependent variables from the previous iteration/time step), it does not make sense to solve eq.
(8.133) to a tight tolerance. Normally, reduction of the sum of absolute residuals for one order
of magnitude suffices.
In order to promote convergence of this solution method, an under-relaxation is often necessary.
This has been introduced in the following manner, by replacing A and b in eq. (8.133) with

1  1 
A A+ D and b b+ D (n-1) , (8.135)
 

respectively, where  is an under-relaxation factor with a value between 0 and 1 (typically

0.9), D is a diagonal matrix consisting of the diagonal elements of matrix A (coefficients aP)
and  (n–1) is the dependent variable vector from the previous iteration/time step. In addition,
this enhances the diagonal dominance of the linearised equations, which improves the rate of
convergence of most iterative linear equation solvers. Also, if the source term is a function of
P, a suitable linearisation may provide another positive contribution to the central coefficient.

The segregated solution strategy employed enables reuse of the same storage for matrix A and
vector b for all dependent variables , thus requiring only 6N storage locations in a 2D case (8N
storage locations in a 3D case).

It is also important to recall that the fully implicit time differencing used avoids stability-related
time step restrictions. In principle, it allows any magnitude of the time step to be used, and in
practice it is limited only by the required temporal accuracy. The number of iterations per time
step depends on the size of the time increment t; for smaller t fewer iterations are required to
obtain the solution at the time tm. If only a steady state is sought, the temporal accuracy is not
important and the solution can be obtained by using a single, infinite time step.

147
 8. FINITE VOLUME DISCRETISATION

In case of fluid flow calculations it may, however, be useful to mention the following:
(1) Due to the presence of the convective transport fluid flow equations are always non-
linear.
(2) Due to the deferred-correction implementation of the discretised convection terms, which
employs the asymmetric upwind scheme, the coefficient matrix of the linearised momentum and
energy equations (8.133) is always asymmetric and the ICCG method cannot be used (one can
employ the Gauss-Siedel method, or some variant of the incomplete Cholesky preconditioned
conjugate gradient method applicable to systems with asymmetric matrices, e.g. CGSTAB).
(3) The under-relaxation factors are typically: for the momentum equations v = 0.7 – 0.8,
for the pressure correction equation p'  1 – v, and for the energy equation T = 0.9 – 1.
(4) From the point of view of a segregated solution procedure the existence of the
continuity equation is very convenient, because it is a single scalar equation which links all
velocity components and 'takes care' of the inter-equation coupling. Consequently, unlike
thermo-elastic stress analysis where the employment of a multigrid acceleration procedure is
crucial for a good convergence rate, coupling among velocity components is not a bottle-neck
for convergence, and the problem of convergence moves towards resolving non-linearities.
Finally, it may be useful to explicitly state the way the pressure-correction equation is
incorporated into the solution algorithm:

Step 1. Provide initial values of dependent variables (at the time to).
Step 2. Assemble and solve equation (8.133) for the velocity components, employing the
currently available pressures and mass fluxes.
Step 3. Assemble and solve equation (8.133) for the pressure correction and use calculated
values to correct mass fluxes, velocity components and pressure [eqs. (8.129) and (8.130)].
Step 4. Assemble and solve equation (8.133) for internal energy and obtain temperature.
Step 5. Return to Step 2 and repeat until the sum of the absolute residuals for all equations
has fallen by a prescribed number (typically three) orders of magnitude.
Step 6. Advance the time by the time increment t and return to Step 2; repeat until the
prescribed number of time steps is completed.

148
 8. FINITE VOLUME DISCRETISATION

8.4 Summary

The main features of the present FV discretisation method are discussed and

• a fully implicit temporal differencing scheme, allowing transient problems to be solved by

using arbitrary large time steps without the numerical stability problems associated with explicit
methods, and

• a second order accurate (on a uniform mesh) fully conservative spatial differencing, enabling
reasonably accurate and physically meaningful results to be obtained on a coarse numerical
mesh are selected.
Then, the 2D momentum and energy balance equations are discretised for cases which can be
'meshed' by a non-uniform Cartesian mesh. This resulted in a 5-diagonal system of 3 x N
generally non-linear algebraic equations for each time step.

The adopted procedure of solving large systems of non-linear algebraic equations incorporates:
• a segregated iterative solution algorithm offering a very low memory requirement which
enables large problems to be solved on small computers; in addition, due to its iterative nature,
it is directly applicable to non-linear solids,

• an efficient conjugate gradient linear equations solver which is considered to be one of the
most powerful methods for the solution of linear algebraic equations, which involves little, or,
sometimes, no memory overheads.

149
 8. FINITE VOLUME DISCRETISATION

8.5 Three-Dimensional Problems in Domains of Arbitrary Shape

In the previous section a 2D FV discretisation involving geometrically simple situations was

described. However, the great majority of engineering problems involve solid bodies of
complex shape. In order to solve such problems, the space discretisation by a structured
Cartesian mesh has to be replaced by a more general one.

In this Section a second-order accurate numerical method will be presented which enables a
solution domain to be discretised by polyhedral cells of a completely arbitrary shape, thus
enabling geometrically complex problems to be solved accurately and efficiently.

Due to the time constraint, the energy and momentum equations for a thermo-elastic solid will
only be considered. The discretisation of the fluid equations normally requires discretisation of
the convective terms [see Sections 8.1.4 and 8.2.2] and pressure terms (employing continuity
equation, Section 8.2.2), and will not be considered here.
In case of thermo-elastic solids, the generic transport equation (7.39) can be simplified as:

t  
B dV =   grad   ndS +  q S ndS +  qV dV , (8.136)
V S S V

where  stands for the displacement u and temperature T, and all other terms are defined in table
1 [Section 7.5].

8.5.1 Numerical method

Before eq. (8.136) is integrated several important decisions have to be made concerning: (a) the
choice of vector and tensor components, (b) space, time and equation discretisation procedure
and (c) variable arrangement. An appropriate decision about those options is decisive for the
stability, conservativeness and efficiency of a numerical method. In the present method the
following choices are made.

a) Vector and tensor components related to a global Cartesian coordinate system are
preferred because they lead to a strong conservation form of momentum equations,
without terms requiring a smooth numerical mesh.
b) The space is discretised by an arbitrary unstructured mesh, because it is by now generally
accepted that only unstructured meshes are capable of describing accurately complicated
3D geometries. The time is subdivided into a number of unequal time intervals. Spatial
gradients are interpolated by a second-order accurate scheme, while for the temporal
discretisation a fully implicit time differencing is employed.
c) All dependent variables are stored at the cell centre, i.e. a collocated variable arrangement
is used, which requires only one set of control volumes to be generated.

150
 8. FINITE VOLUME DISCRETISATION

Details about the implementation of the above options into the present method will now be
explained.

Space and time discretisation

In order to obtain discrete counterparts of the eq. (8.136), the solution domain is discretised into
a finite number of contiguous control volumes (CV) or cells and time is subdivided into an
arbitrary number of time steps tm. The computational points (nodes) are placed at the centre of
each CV, while boundary nodes, needed for the specification of boundary conditions, reside at
the centre of boundary cell faces. The control volume is defined by the coordinates of its
vertices and can be of an arbitrary polyhedral shape, i.e. it can have an arbitrary number of cell
faces Sj (j = 4,5,6,...) (Fig. 8.29). In the same problem control volumes of different topologies
can be used.
N1
N2
n1

n2 S1
S2

P
nn Sn
Nn
Sj

z nj
xP
dj
Nj
y
x

Figure 8.29 Control volume of an arbitrary polyhedral shape.

After that eq. (8.136) for each control volume is written as follows
n n

B dV =    grad   ndS +   qS  ndS +  qV dV ,
t  
(8.137)
V j =1 Sj j =1 Sj V

Rate of change Diffusion Source

where n is the number of cells which share cell faces with cell P (so called nearest neighbours)
and  and B in the case of momentum equations are now related to the Cartesian displacement
components ui (i = 1,2,3).

151
 8. FINITE VOLUME DISCRETISATION

Eq. (8.137) has three distinct parts: transient rate of change, diffusion and source. This equation
is exact, i.e. no approximation has been introduced so far. However, in order to evaluate
integrals in the above equation, a way of calculating integrals and a distribution of dependent
variable  and physical properties in space and time have to be assumed.

Calculation of integrals
Both volume and surface integrals are calculated by using the mid-point formula.

Spatial distribution

Here, the following linear spatial distribution is employed (although, in principal, any other
distribution can be used)
 (x,t ) =  P ( t) + gP ( t)( x  x P ) (8.138)
where  stands for dependent variable , physical properties of continuum or gradient of , xP
is the position vector of point P and gP(t) is an approximation of gradient of  at point P.

The unknown vector gP is calculated by ensuring a fit to a set of sampling points usually
consisting of the nearest neighbours of point P, i.e. by solving the following set of equations
d j g P =  N j   P ( j = 1,...,n) (8.139)

where d j = x N j  x P is the distance vector joining point P with its neighbour Nj (Fig. 8.29). To
solve this over determined set of equations the least square method is used, resulting in the
following normal equation
-1
DgP = h i.e. gP = D h (8.140)

where

 dTj ( N )
n n
D=  dTj  d j and h= j
P (8.141)
j =1 j =1

It is important to notice that for all 's and for assumed linear distribution, D is a symmetric 3x3
matrix, whose coefficients depend solely on the mesh geometrical properties and hence, for a
given cell, a single evaluation of D1 suffices.
Formula (8.138) is used to calculate values of  at the faces j, necessary for evaluation of the
surface integrals in (8.137), leading to a second-order symmetric formula

j =
1
2 ( 1
) ( ) ( )

 P +  N j +  (grad ) P x j  x P + (grad ) N x j  x N j
2 j 
(8.142)

where xj is the position vector of the cell-face centre (Fig. 8.29).

152
 8. FINITE VOLUME DISCRETISATION

Temporal distribution
As far as temporal distribution of is concerned, according to the adopted fully implicit time
discretisation all values of dependent variables are expressed as

( x,t) =
(x,tm ) =
m (8.143)

where m is the time step counter and tm is the current instant of time, except in the transient
term where a linear variation of in time is assumed. Henceforth the time step counter will
normally be omitted, and all values of will refer to the current time tm, unless indicated
otherwise.

Rate of change

The transient term in eq. (8.137) is approximated as


t  B dV 
1
t m
(

) ( m1

 BV P  B V P  ) (8.144)

V

where the meaning of B and B m1 for different is given in Table 8.1.

m1
Table 8.1: The meaning of B and B in eq. (8.144).
m1
B B
m1 m1
T cT c T
ui
1
t m ( m1
ui  ui ) 1
t m1 ( m1
ui  ui
m2
)
Diffusion flux
Diffusion flux Dj of through the cell face j can be approximated by

Dj =   
grad
 ndS  
j (grad
) j n j S j ,
*
(8.145)
Sj

where  j stands for the cell-face mean value of the diffusivity, obtained by using second-order
formula (8.142), and the value of (grad ) j is constructed in the following manner
*

N  P grad  d 
(grad )*j = (grad ) j +  j

j
n , (8.146)
 dj dj  j


where the overbar means an arithmetic average between nodes P and Nj straddling the cell face j
and value of (grad ) j is calculated by using second-order formula (8.142).
The gradient based on expression (8.142) is second order space-centred and as such cannot
sense the oscillations which have the period twice the characteristic length of the numerical
mesh (see expression 8.125 and Fig. 8.27 in Section 8.2.2). As a result, once induced, the
unphysical oscillatory profile remains superimposed onto the otherwise smooth spatial variation

153
 8. FINITE VOLUME DISCRETISATION

of the dependent variable. Because of this, a third-order diffusion term (term in [] brackets) is
introduced into expression (8.146) for the cell-face gradient of variable . This additional
diffusion transport smoothes out these unphysical oscillations and, once a converged solution is
obtained, it will vanish if the spatial variation of is linear or quadratic. Otherwise, its
magnitude is proportional to the truncation error of the second-order scheme, so it does not
affect the second-order behaviour of presented discretisation practice.

The contribution to the diffusion flux coming from the second term in eq. (8.146) is treated
implicitly. The rest of the diffusion flux is due to component of the gradient normal to the
distance vector dj (the 'cross-diffusion') and it vanishes when the grid is orthogonal, and is small
compared to the so called 'normal diffusion' if the grid non-orthogonality is not severe. For this
reason it is treated explicitly. This simplifies the coefficient matrix of the resulting set of linear
algebraic equations since only the contributions of the nearest neighbours of node P are then
treated implicitly.

Source term
The part of the source term coming from the true, volume sources is integrated by assuming a
linear variation of the source over the control volume
QV =  qV dV  qV VP ,
P
( ) (8.147)
V

while the part coming from the stress tensor in momentum equations is discretised in the same
manner as diffusion term

 [μ(gradu) ]
n n

 + [divu  (3 + 2μ) (T  T0 )]I  ii  n j S j (8.148)
T
Qu iS = qS  ndS 
j =1 S j j =1
j 

where ii (i = 1,2,3) are the Cartesian base vectors. Note that qTS = 0.
It is often possible to linearise volumetric source term Q V and add a part of it to the central
coefficient of the resulting coefficient matrix which will then increase its diagonal dominance,
thus enhancing the convergence of the resulting linear algebraic system. It may also be
beneficial, from a stability and efficiency point of view, to take some parts of the surface source
Q S into account implicitly, as was done in Section 8.2.

Resulting algebraic equations and solution procedure

After assembling all the terms featuring in eq. (8.137), it can be written in the form of the
following non-linear algebraic equation which links the value of dependent variable at the
control volume centre with its values at the points in the neighbourhood
n
aP P   aj Nj = b (8.149)
j =1

154
 8. FINITE VOLUME DISCRETISATION

where
Sj
a
j = 
j ,
dj
n E
VP
a
P =  a
j + tm
, (8.150)
j =1
n  grad
d j 
m1 m1
E
VP
b
=  
j (grad
) j  n j 
S + Q
V + Q
S +
j =1
dj j tm

m1
and E and E are defined in Table 8.2.

m1
Table 8.2: The meaning of E  and E in eq. (8.150).
m1
E E
T c c m1T m1
ui

t m
 m1
t m
ui +

u
tm1 i (
m1 m2
 ui )
Eq. (8.149) is of the same form as eqs. (8.101), (8.103) and (8.108), and it is easy to verify that
the coefficients (8.150) reduce to the corresponding ones from Section 8.2 [(8.100), (8.102) and
(8.107)] if a Cartesian space discretisation is employed. The type and implementation of the
initial and boundary conditions is analogous to those described in Section 8.2, while the
segregated solution procedure is entirely applicable to the solution of the system of equations
(8.149) with coefficients and sources (8.150).

8.5.2 Conclusions
The numerical method presented here can be applied to a wide range of problems in continuum
mechanics. Although very general as far as its applicability to different geometrical and
physical situations is concerned, the method is extremely simple. Namely, at first sight an
arbitrary CV that can be used for space discretisation (Fig. 8.29) might suggest that a method
which can use such discretisation must involve very complicated mathematics. The fact is,
however, that all one has to deal with are vectors and tensors in a global Cartesian coordinate
frame (unless someone wishes to do otherwise). As far as physics is concerned, the method
reduces to a straight forward application of the first principles to a control volume.

The use of cells of different topologies that can be combined in an arbitrary way is the most
general approach to space discretisation and greatly facilitates the mesh generation for complex
geometrical configurations when an accurate description of boundaries is desired. This feature
is particularly useful for adaptive mesh refinement when injection and depletion of
computational points is straightforward and at the same time should help automatic mesh
generation, which is rapidly becoming a very important issue in computational continuum
155
 8. FINITE VOLUME DISCRETISATION

mechanics. The discretisation presented facilitates the use of such meshes. It is performed in
real space, using global Cartesian coordinates, thus eliminating the need for coordinate
transformation, which makes it easy to understand and implement. The discretisation promotes
the conservation of momentum and energy on a discrete level, which, combined with a fully
implicit time differencing, enables physically meaningful results to be obtained on very coarse
meshes and for large time steps.

156
 9. INTRODUCTION TO THE FINITE ELEMENT METHOD

Chapter 9

An Introduction to the Finite Element Method

9.1 Background
As discussed in chapters 1–8, continuum mechanics involves the solution of equations (e.g.
conservation of momentum) defined over a body, subject to initial and boundary conditions.
When the geometry of the body of interest is simple, analytical methods can be used (e.g.
beam theory) and exact solutions may often be obtained. Such methods cannot easily
be applied to bodies of complex shape. However, by dividing the body into small, finite
‘elements’ of regular shape and applying the laws of continuum mechanics to the collection
of elements, an approximate solution to the problem can be obtained. As the number of
elements used is increased (and their size decreased) we approach the exact continuum
mechanics (infinitely small elements) solution.
The development of the finite element method is often attributed to Clough and co-
workers in the 1960s at Boeing in the US, who coined the term finite element, and used
the method to analyse aircraft structures. Argyris at Imperial College was also one of the
pioneers in this field. (A very brief and personal history of the finite element method is
provided by Hughes of Stanford University at www.usacm.org/T.J.Hughesacceptanceremarks.htm).
The solution of a finite element problem can involve many thousands of repetitive
calculations and therefore most problems must be solved by computer. With the increasing
availability and affordability of computers in the 1980s and 1990s, the finite element method
expanded rapidly. Nowadays, most structural design makes use of the finite element method,
as it can significantly decrease the time and cost to market of a product, by reducing the need
for expensive prototype testing. Many commercial finite element codes exist, e.g. ABAQUS
(www.abaqus.com), ANSYS (www.ansys.com), ADINA (www.adina.com) and DYNA (www.lstc.com),
most of which are available within the mechanical engineering department. The largest
finite element problems can involve millions of variables (degrees of freedom) and can take
weeks to solve even on the fastest computers currently available.
The application of the finite element method was originally limited to static solid me-
chanics problems but has since been extended to the study of many field problems, e.g.
heat transfer, fluid mechanics and impact and dynamics problems. In this part of the course
we focus on the solution of solid mechanics problems, but the methods presented can be
relatively easily extended to deal with more complicated problems.
Some of the content of these notes is based on the textbook by Belytschko, Liu and
Moran: Nonlinear Finite Elements for Continua and Structures, published by Wiley and
sons, 2000. Another useful reference book, which is perhaps easier to follow for the beginner
but uses somewhat different notation from these notes, is The Finite Element Method, by

157
 9. INTRODUCTION TO THE FINITE ELEMENT METHOD

Zienkiewicz and Taylor, 2000. This book is in three volumes, the most relevant to this
course being Volume 1: The Basis.

9.2 Description of the Boundary Value Problem

In what follows we assume small displacements and rotations so we do not distinguish

between a deformed and undeformed configuration and the stress matrix is unique and
symmetric. However, there is no major difficulty in generalising the method to account for
large displacements and rotations.
fsn
Γt
ub
V Γu

Figure 9.1, Body V subject to applied displacements and tractions

The finite element method is a displacement-based method in that it involves the determi-
nation of the unknown displacement field, u(x). We consider a body, V , with boundary
Γ ≡ Γu + Γt as shown in Fig. 9.1. Displacements, ubi , are specified on the displacement
boundary, Γu , and tractions, fisn , on the traction boundary.
The key equations are the conservation of linear momentum equations and the traction
boundary conditions,
σji,j + ρfib − ρüi = 0, (9.1a)
σji nj = fisn on Γt . (9.1b)
We seek the displacement field, u, which satisfies these equations as well as the displacement
boundary conditions, u = ub on Γu . Equations 9(a) and (b) are sometimes known as the
‘strong’ form of the momentum balance. We shall convert this strong form to the ‘weak
form’, also known as the principle of virtual work.

9.3 Principle of Virtual Work (See also Section 5.4)

Consider a virtual displacement field, δui . This field may take any value provided it is
continuous and bounded (not infinite) over V and δui = 0 on the displacement boundary,
Γu (see section 5.4.1). Multiply the linear momentum equation, Eq. 9.1(a) by this virtual
field, δui , and integrate over the body V to get:

δui (σji,j + ρfib − ρüi )dV = 0. (9.2)
V

158
 9. INTRODUCTION TO THE FINITE ELEMENT METHOD

Consider the first term in the above equation:

 
∂
(σji,j δui )dV = (σji )δui dV
V V ∂xj
   
∂ ∂
= (σji δui )dV − σji δui dV
V ∂xj V ∂xj

Using Gauss’s theorem (Eq. 2.80) we can rewrite the first term on the RHS of the above
equation as follows:  
∂
(σji δui )dV =  · (σδu)dV
V ∂xj V

= (σδu)ndΓ
Γ

where n is the unit normal to the surface Γ. Then

 
(σδu)ndΓ = σji δui nj dΓ
Γ Γ

= σji δui nj dΓ,
Γt

since δui = 0 on Γu . Therefore the first term in Eq. 9.2 becomes:

  
(σji,j δui )dV = σji δui nj dΓ − σji δui,j dV
V Γt V
  (9.3)
sn
= δui fi dΓ − σji δui,j dV .
Γt V

Insert Eq. 9.3 into Eq. 9.2 and multiply by −1 to get,

 
b
(σji δui,j − δui ρfi + δui ρüi )dV − δui fisn dΓ = 0. (9.4)
V Γt

Equation 9.4 is known as the weak form of the momentum balance equation or the principle
of virtual work. In Eq. 9.4 δui is the virtual displacement field which may take any value
provided it is zero on Γu .
We have shown that the strong form of the momentum equation, Eq. 9.1 leads to the
weak form, Eq. 9.4. It can also be shown that the inverse applies—if Eq. 9.4 is satisfied for
all virtual displacement fields, δui , then the displacement field ui satisfies Eq. 9.1(a) and
the corresponding stress fields satisfies the traction boundary conditions, Eq. 9.1(b). Note
that the unknown in Eq. 9.4 is the displacement field, ui (σij depends on ui through the
strain εij which is the derivative of ui .) It will be seen that the virtual field δui does not
appear in the final form of the finite element equations.

159
 9. INTRODUCTION TO THE FINITE ELEMENT METHOD

9.4 Finite Element Discretisation

We first simplify the weak form equation by assuming equilibrium, üi = 0, taking body
forces as zero, fib = 0, and invoking symmetry of the stress tensor, σji = σij . The equation
becomes:  
σij δui,j dV − δui fisn dΓ = 0. (9.5)
V Γt

Figure 9.2, Body V divided into finite elements

To solve Eq. 9.5 we discretise the body V , as shown in Fig. 9.2, i.e. we divide it up
into small discrete finite sized elements (as opposed to the continuum approach where the
material is continuous and material ‘elements’ are infinitesimally small). Along the sides
of each element we designate nodes which define the shape of the element (see Fig. 9.2).
The elements in Fig. 9.2 are shown as squares (apart from the surface elements which are
triangles or quadrilaterals). In general the elements may take any shape, though as will
be seen, for most finite element formulations the accuracy of the method is highest when
the elements are squares (or rectangles). The full collection of nodes and elements which
represent the body, V , is known as the finite element mesh.
We wish to obtain the unknown displacement field, ui (x), within the domain V . We
define the displacements using shape functions (also known as interpolation functions) as
follows:
Nn

ui (x) = Na (x)Uia , (9.6)
a=1

where Uia are the discrete (unknown) nodal displacements, Na (x) are the (known) shape
functions, which we will discuss further later, and Nn is the total number of nodes in the
model.
The virtual displacements δui are interpolated using the same shape functions:

Nn

δui (x) = Na (x)δUia , (9.7)
a=1

where δUia are the virtual nodal displacements.

160
 9. INTRODUCTION TO THE FINITE ELEMENT METHOD

In the weak form we require derivatives of displacements and virtual displacements.

Again the shape functions are used:
Nn

δui (x) = Na (x)δUia . (9.8)
a=1

Nn

∂ ∂
⇒ δui (x) = Na (x)δUia
∂xj x
a=1 j
(9.9)
Nn

= Na,j (x)δUia
a=1

The virtual work equation, Eq. 9.5, then becomes:

 N   N 
 n n

σij Na,j δUia dV − fisn Na δUia dΓ = 0

V a=1 Γt a=1

Nn
   
⇒ δUia σij Na,j dV − Na fisn dΓ =0 (9.10)
a=1 V Γt

The above equation must hold for any virtual displacements δUia . Therefore, we can write,
 
σij Na,j dV − Na fisn dΓ = 0
V Γt

or  
σij Na,j dV = Na fisn dΓ (9.11)
V Γt

In Eq. 9.11, i ranges from 1 to N DF , where N DF is the number of spatial degrees of

freedom (usually two or three) and a ranges from 1 to Nn . Therefore Eq. 9.11 is a vector
equation with N DF × Nn terms. The term on the right hand side of Eq. 9.11 is often
referred to as the external nodal force vector (it is due to the externally applied tractions)
and the term on the left hand size is referred to as the internal nodal force vector. The
finite element equation then simply becomes,

internal nodal forces = external nodal forces. (9.12)

9.5 Linear elastic finite element model

The finite element equation, Eq. 9.11, may be further simplified by using the constitu-
tive law to write the stress in terms of the unknown nodal displacements. We specialise to
the case of linear elasticity, though there is not too much difficulty in dealing with non-linear

161
 9. INTRODUCTION TO THE FINITE ELEMENT METHOD

material behaviour. (In the latter case the finite element equations become non-linear and
cannot be solved in a single step (iteration)).
For linear elasticity, following Eq. 6.12,

σij = cijkl εkl (9.13)

and using the small strain-displacement relationship, Eq. 4.35, we get,

1
σij = cijkl (uk,l + ul,k )
2
= cijkl uk,l (due to symmetry of cijkl ).
Now
Nn

uk = Nb Ukb ,
b=1

Nn

⇒ uk,l = Nb,l Ukb ,
b=1
so
Nn

σij = cijkl Nb,l Ukb . (9.14)
b=1

Therefore we can rewrite Eq. 9.11 as,

 
Nn 
cijkl Nb,l Ukb Na,j dV = Na fisn dΓ.
V b=1 Γt

Rewrite as a matrix equation

  
cijkl Nb,l Na,j dV Ukb = Na fisn dΓ (9.15)
V Γt

or
Kiakb Ukb = fia
or
KU = f , (9.16)
where U is the vector of unknown nodal displacements (e.g. U21 is the displacement at node
1 in the 2 (y) direction), f is the nodal force vector and K is the finite element ‘stiffness’
matrix (note the analogy with the equation for a spring, F = kx, where F is the force in
the spring, k the spring stiffness, and x the displacement of the spring).
The solution to Eq. 9.16 is obtained by any appropriate numerical method, which
essentially involves the inversion of the stiffness matrix, i.e.

U = K−1 f .

162
 9. INTRODUCTION TO THE FINITE ELEMENT METHOD

Equation 9.16 is the finite element equation for a linear analysis—linear geometry (small
displacements and rotations) and linear material response. The formulation for non-linear
response is similar. It will take the same form as Eq. 9.16 but the stiffness matrix will not
be constant (as in the linear case) but will depend on the applied load or displacement.
The incorporation of dynamic effects (via the ρüi term in the momentum equation) is also
relatively straightforward.

9.6 The Element and Global Stiffness Matrices

The main task in a finite element analysis is the calculation and inversion of the finite
element stiffness matrix. In a non-linear analysis this operation may be carried out many
thousands of times. The matrix K in equation 9.16 is often referred to as the ‘global’
stiffness matrix. The term Kiakb can be interpreted as the stiffness at node a in the degree
of freedom i due to movement of node b in degree of freedom k.
Each element in the finite element mesh contributes an amount to the global stiffness
matrix and the amount it contributes is called the ‘local’ stiffness matrix. The finite element
procedure operates on the principle that all calculations carried out at an element level require
only information from the element itself, which considerably streamlines the computational
analysis as each element can be treated independently and the calculations carried out on
each element are identical. This feature is achieved through the use of the shape functions,
Na , whereby the value of the variables within an element depend only on the values at the
nodes of that element.

9.6.1 Definition of finite element shape functions

The global shape functions Na (x), used to interpolate the displacement field, are
constructed by adding together the local (element) shape functions, Nae (x). For now we
focus on the latter. The requirements for the element shape functions Nae are that
(i) Nae = 1 at node a
= 0 at all other nodes
N
 EL
(ii) Nae = 1 where N EL is the number of nodes per element
a=1
(iii) Nae = 0 outside element e.
Within an element the displacement u is given by

N
 EL
ui (x) = Nae (x)Uia . (9.17)
a=1

For the rest of this sub-section the superscript e is dropped for convenience.

163
 9. INTRODUCTION TO THE FINITE ELEMENT METHOD

Consider an example of a 1-D element with 2 nodes (see Fig. 9.3). The nodes are
designated to lie at x = −1 and x = +1 respectively. (This is known as the parent element.
When carrying out a finite element analysis, the actual elements are mapped to this parent
element (more on this later).)

U1 u(x)

U2

node 1, x = −1 node 2, x = +1

Figure 9.3, 1D linear shape function

The lowest order polynomial that can interpolate a function over the element is a linear
function, as illustrated in Fig. 9.3. In general, for a two noded element, we write
2

u(x) = Na (x)U a
a=1

and the linear shape functions are given by

1
N1 (x) = (1 − x)
2
1
N2 (x) = (1 + x).
2
It is clearly seen that requirement (i) above is satisfied for both shape functions, since
N1 (−1) = 1, N1 (+1) = 0, N2 (−1) = 0, N2 (+1) = 1 and clearly N1 + N2 = 1.
The use of the shape functions in this way is simply a convenient way of defining a
linear interpolation between node 1 and node 2. To find the value of the variable u at any
point within the element,
2
 1 1
u(x) = Na (x)U a = (1 − x)U 1 + (1 + x)U 2 ,
a=1
2 2

e.g.
U1 + U2
u(0) = (as expected).
2
Quadratic interpolation over an element can also be used. For this three nodes are
required, with the third node usually lying halfway between the other two (see Fig. 9.4).
This provides more flexibility in describing the variation of u within an element but also
increases the complexity of the problem as an extra node has been introduced. Note that

164
 9. INTRODUCTION TO THE FINITE ELEMENT METHOD

the 3 node quadratic element in Fig. 9.4 could be replaced by two 2 noded linear elements.
In general quadratic elements are preferred, though there are situations when the use of
quadratic elements can cause problems. For example, in contact problems linear elements
are recommended.
U1

U3 U2

node 1 node 3 node 2

Figure 9.4, 1D quadratic shape function

For the 1-D quadratic element illustrated in Fig. 9.4,

x(x − 1)
N1 (x) =
2
x(x + 1)
N2 (x) =
2
N3 (x) = 1 − x2

For the more common two dimensional elements we can again describe linear and
quadratic shape functions, though in 2D they are actually bi-linear and bi-quadratic (e.g.
see Fig. 9.5). For example, the linear shape functions are given by,

1
N1 (x, y) = (1 − x)(1 − y)
4
1
N2 (x, y) = (1 + x)(1 − y)
4
1
N3 (x, y) = (1 + x)(1 + y)
4
1
N4 (x, y) = (1 − x)(1 + y)
4

Figure 9.5, 2D linear shape function

165
 9. INTRODUCTION TO THE FINITE ELEMENT METHOD

Shape functions of higher than quadratic order are rarely used except for specialised appli-
cations.

9.6.2 Global Stiffness Matrix

We have seen in section 9.5 that the global stiffness matrix may be written as

Kiakb = cijkl Na,j Nb,l dV .
V

The terms Na and Nb in the above are the global shape functions. The global shape
functions are simply the sum over the body of the local shape functions, i.e.
Ne

Na (x) = Nae (x),
e=1

where Ne is the total number of elements in the model.

By rule (iii) the local (element) shape functions are restricted to be zero outside the
element, i.e. Nae (x) = 0 outside the element e. This ensures that the displacements are
continuous throughout the body and that the displacement within an element depends only
on the nodal displacements for that element.
The construction of the global shape functions is best illustrated by a diagram as in
Fig. 9.6. In Fig. 9.6, N2 (x) = N2e=1 (x) + N1e=2 (x).

N2 ( x ) N1e ( x )

1 1 2 2 3 3 4
1 2

Figure 9.6, Global and element shape functions

Note the distinction between local element numbering (shown underlined in Fig. 9.6) and
global element numbering. When setting up a finite element mesh it is the global numbering
that is specified by the user. The finite element program must then keep track of the local
numbering (for a 2 noded element, the local node number is always either 1 or 2).
The global stiffness matrix is written using the local shape functions as
 Ne
 Ne

e f
Kiakb = cijkl Na,j Nb,l dV
V e=1 f =1
(9.18)
Ne 

e e
= cijkl Na,j Nb,l dV .
e=1 Ve

166
 9. INTRODUCTION TO THE FINITE ELEMENT METHOD

1 2
The cross terms (e.g. N2,j × N1,l ) disappear because of the property that Nae (x) = 0 and
Naf (x) = 0 outside the element e or f respectively.
We can then rewrite the stiffness matrix as
Ne

e
Kiakb = kiakb ,
e=1

e
where kiakb is the element stiffness,

kiakb = cijkl Na,j Nb,l dV , (9.19)
V

(dropping the superscript e for convenience).

As discussed previously, one advantage of the above form is that, although the global
stiffness matrix is made up of contributions from the element stiffness matrix, the element
stiffness matrix (sometimes called the local stiffness matrix) depends only on quantities
within that element, which makes programming of the method relatively simple. Construc-
tion of the global stiffness matrix involves evaluation of Eq. 9.19 for each element in the
mesh, followed by the addition of each local stiffness into the global stiffness matrix (usually
called the ‘assembly’ process). The assembly process will be discussed in more detail in
subsequent sections.

9.6.3 Computation of Element Stiffness Matrix

9.6.3.1 Use of matrix (Voigt) notation

So far we have used tensorial notation to denote the stiffness matrix and the finite
element equation. In some cases the equations are more easily represented if we use the
equivalent ‘matrix’ notation (sometimes referred to as Voigt notation).

9.6.3.1.1 Stress strain relation

Consider the constitutive law for an isotropic linear elastic material (Eqs. 6.12 and
6.15),

σij = cijkl εij = λεkk δij + 2μεij (9.20)

We can write this in matrix form as follows. The stress tensor is represented as a vector,
⎡ ⎤
σ11
⎢ σ22 ⎥
⎢ ⎥
⎢ σ33 ⎥
σ=⎢ ⎥ (9.21)
⎢ σ12 ⎥
⎣ ⎦
σ13
σ23

167
 9. INTRODUCTION TO THE FINITE ELEMENT METHOD

and then the constitutive equation may be rewritten as

⎡ ⎤ ⎡ ⎤⎡ ⎤
σ11 λ + 2μ λ λ 0 0 0 ε11
⎢ σ22 ⎥ ⎢ λ λ + 2μ λ 0 0 0 ⎥ ⎢ ε22 ⎥
⎢ ⎥ ⎢ ⎥⎢ ⎥
⎢ σ33 ⎥ ⎢ λ λ λ + 2μ 0 0 0 ⎥ ⎢ ε33 ⎥
⎢ ⎥=⎢ ⎥⎢ ⎥. (9.22)
⎢ σ12 ⎥ ⎢ 0 0 0 μ 0 0 ⎥ ⎢ 2ε12 ⎥
⎣ ⎦ ⎣ ⎦⎣ ⎦
σ13 0 0 0 0 μ 0 2ε13
σ23 0 0 0 0 0 μ 2ε23

In matrix (Voigt) form this is

σm = Dmn εn , (9.23)

with m, n ranging from 1 to 6. Note the factor of 2 in the shear strain terms in Eq. 9.22
(not all finite element codes follow this convention). This is for convenience and to allow
us to write for example,
σ · ε = σn εn .

9.6.3.1.2 Strain-displacement relation using matrix notation

The strain within an element is determined by differentiating the displacements,

N EL
1 1  
εij = (ui,j + uj,i ) = Na,j Uia + Na,i Uja . (9.24)
2 2 a=1

For example,

N EL N EL
1  a a

ε11 = (Na,1 U1 + Na,1 U1 ) = Na,1 U1a .
2 a=1 a=1

N EL
1 
ε12 = (Na,2 U1a + Na,1 U2a ).
2 a=1

We can write this in matrix form as follows:

⎡ ⎤ ⎧⎡ ⎤ ⎫
ε11 ⎪ Na,1 0 0 ⎪
⎪
⎪ ⎪
⎢ ε22 ⎥ ⎪ ⎢ 0 N 0 ⎥ ⎡ a ⎤⎪ ⎪
EL ⎪ ⎥ U1 ⎪
a,2
⎢ ⎥ N ⎨⎢ ⎬
⎢ ε33 ⎥ ⎢ 0 0 Na,3 ⎥ ⎣ a ⎦
⎢ ⎥= ⎢ ⎥ U2 . (9.25)
⎢ 2ε12 ⎥ ⎪⎢ Na,2 Na,1 0 ⎥ ⎪
⎣ ⎦ a=1 ⎪⎪ ⎣ ⎦ U3a ⎪⎪
2ε13 ⎪
⎪ ⎪
⎪
⎩ Na,3 0 Na,1 ⎭
2ε23 0 Na,3 Na,2

We can expand the sum in the strain-displacement relation above to get (for a 4 noded
element):

168
 9. INTRODUCTION TO THE FINITE ELEMENT METHOD

Equation 9.26 is for a 3D analysis. It can be simplified if we consider 2D approximations,

(e.g. for plane strain ε13 = ε23 = ε33 = 0, U3a = 0).
In general the strain-displacement equation of Eq. 9.25 can be written as:

εm = Bmp Up , (9.27),

where B is the matrix of shape function derivatives, as illustrated in Eq. 9.26. In Eq. 9.27
m ranges from 1 to N ST R, where N ST R is the number of strain components (6 in Eq.
9.26) and p ranging from 1 to N EL × N DF , where N EL is the number of nodes per
element (4 in Eq. 9.26) and N DF is the number of degrees of freedom per node (3 in Eq.
9.26).

9.6.3.1.3 Element stiffness matrix using matrix notation

If we recall the original form of the weak form expression (Eq. 9.5), the contribution
of the stiffness term for a linear elastic material may be written as
 
σij δui,j dV = cijkl uk,l δui,j dV .
V V

Looking at the integrand of the term on the RHS of the above equation, we can write,

cijkl δui,j uk,l = cijkl δεij εkl .

In ‘matrix’ form,
cijkl δεij εkl = Dmn δεm εn ,

with m and n ranging from 1 to 6 in the general case.

Since
εm = Bmp Up ; δεm = Bmp δUp

169
 9. INTRODUCTION TO THE FINITE ELEMENT METHOD

we can rewrite the above as

cijkl δεij εkl = Dmn Bmp δUp Bnq Uq

(9.28)
= (Bmp Dmn Bnq )δUp Uq .
By consideration of Eq. 9.10 and 9.15, it may be seen that the term in brackets on the
right hand side of Eq. 9.28 may be identified as the global stiffness matrix K. It may further

be seen from Eq. 9.28 that K = V BT DBdV . The element stiffness matrix is then,

k= BT DBdV , (9.29)
V
where D is the elasticity matrix which contains the material properties and B is the matrix
containing the derivatives of the shape functions, (e.g. Eq. 9.26). This is the most
commonly used notation for the finite element stiffness matrix.
To evaluate kiakb we integrate over the element as shown in the above equation. The
value of the shape functions and their derivatives on the parent element are known (see
section 9.6.1) so B is known. The matrix D is made up of material constants (see Eq.
9.22) and so is also known. Integration of the stiffness matrix over the parent element
can therefore be carried out directly. However, as will be discussed, for arbitrarily shaped
elements, a mapping between the parent and the actual element is required, which means
that the stiffness matrix must be integrated numerically. The approach used in the finite
element method is called numerical quadrature and is analogous to the use of the midpoint
rule in finite volume formulations.

9.7 Numerical Quadrature

Consider the one dimensional integral
 +1
I= f (x)dx.
−1
We can replace this integral by a sum,

 +1 N
 QP
I= f (x)dx ≈ wp f (ξp ), (9.30)
−1 p=1

where ξp give the position of the quadrature points (sometimes called gauss points, as the
method is also known as Gaussian quadrature), wp are the ‘weights’ assigned to each gauss
point and N QP is the order of quadrature. If N QP = 1 we have linear quadrature; if
N QP = 2 we have quadratic quadrature etc.
If we have N QP quadrature points it may be shown that we can integrate exactly a
polynomial of order 2 × N QP − 1, i.e. with 1 gauss point we can integrate a linear function
exactly and if N QP = 2 we can integrate a cubic function.

170
 9. INTRODUCTION TO THE FINITE ELEMENT METHOD

The coordinates and weights of the gauss points corresponding to a particular order
of integration are given in any standard text book on numerical methods. They may be
evaluated by integrating simple expressions whose results are known and comparing with the
result obtained from Gauss quadrature. This is illustrated below for the case of N QP = 1.
9.7.1 First Order Quadrature
For N QP = 1 we consider integrating a constant function f (x) = 1 and a linear function
f (x) = x over the ‘parent’ element.
 1 1

1dx = 2; wp f (ξp ) = w1 × 1 = w1
−1 p=1

⇒ w1 = 2

 1 1

xdx = 0; wp f (ξp ) = w1 × ξ1 = 2ξ1
−1 p=1

⇒ ξ1 = 0.

Therefore the general expression for one gauss point N QP = 1 is

 1
f (x)dx = 2f (0), (9.31)
−1

i.e. the midpoint rule. Equation 9.31 will exactly integrate functions of the form

f (x) = a0 + a1 x.

9.7.2 Second and Higher Order Quadrature

To integrate higher order polynomials we need higher order integration. For N QP = 2
√ √
it can be shown that w1 = w2 = 1 and ξ1 = −1/ 3, ξ2 = 1/ 3. Table 9.1 shows the
location of the gauss points and the magnitudes of the weights for values of N QP from 1
to 4.

NQP ξ w
1 0 2
2 −1/√3 +1/√3 1 1
3 − 6 / 10 0 − 6 / 10 5/9 8/9 5/9
4
− (3 + 4.8 ) / 7 − (3 − 4.8 ) / 7 (3 − 4.8 ) / 7 (3 + 4.8 ) / 7 1/2 − √30/36 1/2 + √30/36 1/2 + √30/36 1/2 − √30/36

Table 9.1, Location of gauss points and magnitudes of weights for N QP = 1 to 4.

171
 9. INTRODUCTION TO THE FINITE ELEMENT METHOD

9.7.3 Integration over a 2 Dimensional Space

For the more general 2D case we are integrating functions of the form,
 1  1
I= f (x, y)dxdy.
−1 −1

We replace the integral I by a sum,

   N QP 
1 1 1 
I= f (x, y)dxdy = wp f (ξp , y) dy
−1 −1 −1 p=1
N
 QP  1 
= wp f (ξp , y)dy
p=1 −1

N QP
N QP 
 
= wp wq f (ξp , ξq ) ,
p=1 q=1

i.e.
 1  1 N
 QP N
 QP
f (x, y)dxdy = wp wq f (ξp , ξq ). (9.32)
−1 −1 p=1 q=1

For NQP = 1, w1 = 2 and ξ1 = 0 (as for 1D) and

I = w1 w1 f (ξ1 , ξ1 ) = 2 × 2f (0, 0) = 4f (0, 0).

For NQP = 2 we have 2 × 2 integration as shown in Fig. 9.7 and Eq. 9.32 is,

2 
 2
I= f (ξp , ξq ),
p=1 q=1

(since wi = 1 for N QP = 2 (see Table 9.1)).

× × × × ×
− − −
× × × ×
− √ √
× × × × ×

− − −

Figure 9.7, Quadrature points for 1 × 1; 2 × 2 and 3 × 3 integration.

172
 9. INTRODUCTION TO THE FINITE ELEMENT METHOD

9.8 Evaluation of Stiffness Matrix

We can now apply Gaussian quadrature to the integration of the element stiffness
matrix. Assuming a finite element mesh with two spatial dimensions,

k= BT DBdV
V
N
 QP N
 QP (9.33)
T
= wp wq B (ξp , ξq )DB(ξp , ξq ).
p=1 q=1

In the above it is assumed that the material properties, D, (λ and μ for isotropic elasticity)
do not vary with position.
We next consider what order of integration is needed to integrate the above equation
exactly. The D matrix is assumed independent of position and does not affect the order of
integration needed. Since the B matrix is made up of the derivatives of shape functions,
the order of B will depend on the order of shape function used, e.g. linear or quadratic.
If we have a linear 2D element the shape functions are bi-linear. Therefore the B matrix
is linear and the integrand of the stiffness matrix, BT DB is quadratic. Therefore we need
N QP = 2. If we are using a quadratic element, then BT DB is quartic (to the power of
4). Therefore in this case, to integrate exactly, we require N QP = 3 (see Fig. 9.7).
The above analysis assumes that the element is equivalent to the parent element. We
will subsequently generalise our analysis to the case of an arbitrarily shaped element.

9.8.1 Output of stress and strain data

To assemble the stiffness matrix using Eq. 9.33 the matrix B is evaluated at the
integration points. Therefore it is relatively simple to evaluate strain at the integration
points using Eq. 9.27 and then stress using Eq. 9.23. Thus finite element codes often
provide as output the stress and strain at the integration points rather than e.g. at the nodal
points (note that the displacement is always obtained at the nodes). However, there is no
difficulty, if required, in extrapolating stress or strain values determined at the integration
points to the nodes or indeed averaging the values over an element, which often provides a
more accurate value for the stress.

9.9 Mapping of Elements to and from the Parent Space

As discussed in Section 9.8 we carry out the integration on the parent element using
numerical quadrature. The next step is to map the result to the actual element (see Fig.
9.8). In what follows the notation s and s1 , s2 will be used to denote coordinates on the
parent element and x and x1 , x2 to denote the mapped element (a 1D or 2D analysis is
assumed throughout though the procedure is identical in 3D).

173
 9. INTRODUCTION TO THE FINITE ELEMENT METHOD

Figure 9.8, Mapping from the parent element to the finite element for 1D and 2D elements.
Consider integration of a function f (x) over a 1D finite element space,

I= f (x)dx
x
By using a change of variables, we can integrate over the parent element, e.g.,
 
dx
f (x)dx = f (x(s)) ds,
x s ds
The term dx/ds is known as the jacobian, J, of the mapping and we can write,
 
f (x)dx = f (s)J(s)ds. (9.34)
X S

To calculate the jacobian we need to know the mapping from s space to x space.
Generally in a finite element analysis, the shape functions are used (i.e. to account for the
element ‘shape’).

(s1, s2)

(x1, x2)

(S 1
2
,S 2
2 )

(X 1
2
, X 22 )
Figure 9.9, Mapping from the s space to the x space.
We write, in the general case, that
N
 EL
x(s) = Na (s)Xia , (9.35)
a=1

174
 9. INTRODUCTION TO THE FINITE ELEMENT METHOD

where Xia indicate the coordinates of the nodes of the mapped element (see Fig. 9.9).
Elements which use the same shape functions to map the nodal coordinates and the
nodal displacements are known as isoparametric elements—the vast majority of elements
do, though not all.
In 2D we define a jacobian matrix, Jij and
N EL  N EL
∂xi ∂  
Jij = = Na (s)Xia = Na,j Xia (9.36)
∂sj ∂sj a=1 a=1

and the jacobian of the transformation from s-space to x-space is the determinant of this
matrix,  
∂xi
J = det .
∂sj
We can then write that
N
 QP N
 QP
k= wp wq BT (ξp , ξq )DB(ξp , ξq )J(ξp , ξq ). (9.37)
p=1 q=1

The value of the shape functions at the mapped gauss points is the same as the value
on the parent element (this is a property of the isoparametric transformation) so they can be
evaluated using the parent element definitions. Note, however, that there are derivatives in
the expression for k. These derivatives are with respect to x, i.e. the mapped coordinates,
while the shape functions are expressed in terms of parent coordinates, s. To carry out the
differentiation we need the inverse of the jacobian matrix, i.e.,
   
∂Na ∂Na ∂si ∂Na −1
= = Jij . (9.38)
∂xj ∂si ∂xj ∂si
It is these mapped derivatives of the shape functions that should be inserted into the
B matrix and used in Eqs. 9.27 and 9.29, if the element of interest is not the same as the
parent element.
9.9.1 Example Calculation of Jacobian for a 4 Noded Element

Consider the 4 noded (linear) element in Fig. 9.10. Writing s1 and s2 as the coordinates
in the parent space, the linear shape functions are given by,
1
N1 (s1 , s2 ) = (1 − s1 )(1 − s2 )
4
1
N2 (s1 , s2 ) = (1 + s1 )(1 − s2 )
4
1
N3 (s1 , s2 ) = (1 + s1 )(1 + s2 )
4
1
N4 (s1 , s2 ) = (1 − s1 )(1 + s2 )
4
175
 9. INTRODUCTION TO THE FINITE ELEMENT METHOD

Figure 9.10, 4 noded parent element in s-space

To integrate over x-space we need the jacobian J of the mapping from s-space to the
x-space (see Fig. 9.11).

Figure 9.11, Mapping from s-space to x space

We write (see Eq. 9.35) that

N
 EL
xi (s) = Na (s)Xia
a=1

where Xia indicate the coordinates of the nodes of the mapped element. (It is easy to
see, using the linear shape functions, that, for example, node 1 on the parent element,
with coordinates, (−1, −1), is mapped to node 1 on the mapped element with coordinates
(X11 , X21 ).)
The jacobian matrix is then
N EL  N EL
∂xi ∂  
Jij = = Na (s)Xia = Na,j Xia
∂sj ∂sj a=1 a=1

The jacobian of the transformation from s-space to x-space is the determinant of this matrix,

J = det (J) .

176
 9. INTRODUCTION TO THE FINITE ELEMENT METHOD

The derivative of the shape functions may be written down immediately, e.g. for the
linear shape functions,

1 1
N1,1 (s1 , s2 ) = − (1 − s2 ) ; N1,2 (s1 , s2 ) = − (1 − s1 )
4 4
1 1
N2,1 (s1 , s2 ) = (1 − s2 ) ; N2,2 (s1 , s2 ) = − (1 + s1 )
4 4
1 1
N3,1 (s1 , s2 ) = (1 + s2 ) ; N3,2 (s1 , s2 ) = (1 + s1 )
4 4
1 1
N4,1 (s1 , s2 ) = − (1 + s2 ) ; N4,2 (s1 , s2 ) = (1 − s1 )
4 4
Consider a regular square element shown in Fig. 9.12.

Figure 9.12, Example used in calculation of jacobian

N
 EL
Jij = Na,j Xia
a=1

Consider each term in turn:

J11 = N1,1 X11 + N2,1 X12 + N3,1 X13 + N4,1 X14

where e.g. X13 is the x1 coordinate of node 3 etc. Then

J11 = N1,1 + 4N2,1 + 4N3,1 + N4,1

1 1
= − (1 − s2 ) + 1 − s2 + 1 + s2 − (1 + s2 ) = 3/2
4 4

Similarly,

J12 = N1,2 X11 + N2,2 X12 + N3,2 X13 + N4,2 X14

= N1,2 + 4N2,2 + 4N3,2 + N4,2
1 1
= − (1 − s1 ) − (1 + s1 ) + (1 + s1 ) (1 − s1 ) = 0
4 4
177
 9. INTRODUCTION TO THE FINITE ELEMENT METHOD

and J21 = 0 and J22 = 3/2. The jacobian matrix is then

 
3/2 0
J=
0 3/2

The jacobian is the determinant of this matrix = 9/4.

Thus in this simple case, the jacobian is simply the ratio of the area of the mapped
element (= 9) to the area of the parent element (= 4). The jacobian is constant over the
element and thus integrating over the mapped element does not require a higher order of
integration than integrating over the parent element, so the integration is exact (see Section
9.7).
−1
Note that unless the element is square J and Jij will not be constant throughout
the element (they will depend on position x). Therefore the order of integration of the
stiffness matrix is higher than that discussed in section 9.7 and integration of the stiffness
matrix will not be exact. Simply having square elements does not of course guarantee
an exact solution to the problem, as the accuracy of the solution is limited by the form
chosen for the shape function—i.e. linear or quadratic variation of displacement. If the
solution to the problem is a linear or quadratic variation in displacement and the geometry
of the body may be represented as a square, then a single square element will provide the
exact solution. However, in general many elements are required in order to construct the
variation in displacement over the whole body by a combination of linear and/or quadratic
distributions over small elements.

9.10 Computation of the Nodal Force Matrix

As discussed in Section 9.5 the global finite element equation is written as (Eq. 9.16),

KU = f

Here f is the nodal force vector, which, as with the stiffness matrix, is made up of contribu-
tions from the individual elements, and in the same way we define a local force vector and
a global force vector.
It is seen from Eq. 9.15 that with applied tractions fisn on Γt the nodal force vector is
given by

fia = Na fisn dΓ. (9.39)
Γt

If body forces exist these are also included as part of the force vector. It may be shown,
starting from the virtual work expression, Eq. 9.4, that the contribution to the force vector
from body forces is given by,

178
 9. INTRODUCTION TO THE FINITE ELEMENT METHOD


fia = Na ρfib dV . (9.40)
V
Typical examples of body forces are weight, where |f | = g and the force acts towards
the ground, and centrifugal force, |f | = ω 2 r, where r is the distance to the axis of rotation
and the force acts along the radius, r. Note that Eq. 9.39 involves integration over the
element surface while Eq. 9.40 involves integration over the element volume. Therefore
the same procedure for integrating the stiffness matrix may be used to integrate the body
force term, but the integration of the traction term, Eq. 9.39, is different. It should also
be pointed out that a different order of integration can be used to integrate the traction
and/or body force term than is used for the element stiffness matrix. For example if the
traction over the element is constant (uniform applied stress) then single point quadrature,
over the element surface is sufficient.
If an applied load may be considered to be a concentrated force, then the force vector
is determined directly, i.e.
fia = Fi . (9.41)
where Fi are the components of the applied force at the node a.

9.10.1 Treatment of Initial Strains

Initial strains (such as thermal strains) may also be included in the nodal force vector. To
see this, it is necessary to examine the original virtual work equation (Eq. 9.11):
 
σij Na,j dV = Na fisn dΓ
V Γt
If initial strains are present, rather than use Eq. 9.13 when substituting for stress in
the above equation, we must use the full form of Hooke’s law:

σij = cijkl (εkl − ε0kl ), (9.42)

where ε0kl are the initial strains (independent of the deformation). A common example of
initial strains are those generated due to a change in temperature. In this case

ε0kl = αΔT δkl ,

where α is the coefficient of thermal expansion, ΔT is the increase in temperature and δkl
is the Kronecker delta.
Inserting Eq. 9.42 into the virtual work equation we get:
 
0
cijkl (εkl − εkl )Na,j dV = Na fisn dΓ
V Γ
  t 
sn
⇒ cijkl εkl Na,j dV = Na fi dΓ + cijkl ε0kl Na,j dV ,
V Γt V

179
 9. INTRODUCTION TO THE FINITE ELEMENT METHOD

i.e., the equivalent nodal force due to the internal strain is given as,

a
fi = cijkl ε0kl Na,j dV .
V

9.11 Assembling the Global Stiffness and Force Matrices

The next stage in the analysis is the assembly of the global system matrices. Each
element will make a contribution to the global stiffness and force matrices. This is best
illustrated by a simple example. Consider a two element truss illustrated in Fig. 9.13. Each
node has two potential degrees of freedom.
F

Figure 9.13, Analysis of a 2 element truss problem

The overall finite element equation for this two element problem is as follows:

⎡ ⎤⎡ 1 ⎤ ⎡ 1 ⎤
K1111 K1121 K1112 K1122 K1113 K1123 U1 f1
⎢ K2111 K2121 K2112 K2122 K2113 K2123 ⎥ ⎢ U21 ⎥ ⎢ f21 ⎥
⎢ ⎥⎢ ⎥ ⎢ ⎥
⎢ K1211 K1221 K1212 K1222 K1213 K1223 ⎥ ⎢ U12 ⎥ ⎢ f12 ⎥
⎢ ⎥⎢ ⎥ = ⎢ ⎥. (9.43)
⎢ K2211 K2221 K2212 K2222 K2213 K2223 ⎥ ⎢ U22 ⎥ ⎢ f22 ⎥
⎣ ⎦⎣ 3 ⎦ ⎣ 3 ⎦
K1311 K1321 K1312 K1322 K1313 K1323 U1 f1
K2311 K2321 K2312 K2322 K2313 K2323 U23 f23

The term Kiakb in the stiffness matrix implies the stiffness at node a in the degree of
freedom i due to movement of node b in degree of freedom k. The stiffness matrix will
be symmetric, Kiakb = Kkaib , and the diagonal terms will be positive or zero (this is true
in general, not just for this example). Since node 2 is common to two elements (see Fig.
9.13) , there will be a contribution from element 1 and element 2 to all terms which have
1 2
a and b = 2, e.g. K1212 = k1212 + k1212 . Note also that since node 3 and node 1 are not
connected, all terms of the form Ki3k1 = Ki1k3 = 0. For this problem assembly of the
global stiffness matrix is relatively trivial, but the principle is the same for all FE meshes.
The assembly of the force vector is also trivial in this case. The only force is a con-
centrated force at node 2 in the 2 direction. Therefore f22 = F and all other nodal forces

180
 9. INTRODUCTION TO THE FINITE ELEMENT METHOD

are zero. There are nodal reaction forces at node 1 and 2, but as will be seen in the next
section, these are not included in the final system of equations.

9.12 Application of Displacement Boundary Conditions and Solution

of the Finite Element Equation
In the problem illustrated in Fig. 9.13, U11 = U21 = U13 = U23 = 0. Therefore we do
not need to solve for these variables and we can remove the columns in the stiffness matrix
corresponding to these degrees of freedom in Eq. 9.43. It can be shown that we can also
remove the rows corresponding to these degree of freedom from the matrix equation and
Eq. 9.43 reduces to   2   
K1212 K1222 U1 0
2 = .
K2212 K2222 U2 F
The problem is then solved by inverting the 2 × 2 matrix and solving for U12 and U22 .
If required the nodal reaction forces, f11 , f21 , f13 , f23 may then be obtained from the full
finite element equation, Eq. 9.43. (The equations corresponding to these unknowns were
removed when the boundary conditions were applied.) Note that if the rows corresponding
to the displacement boundary conditions are not removed the stiffness matrix cannot be
inverted (it is singular). Physically, this corresponds to the situation when a load is applied
to an unsecured body which can give rise to an infinite displacement.

9.13 Steps in a static linear finite element analysis

1. Define area of interest—number of spatial dimensions, geometrical boundaries of the
body, number of degrees of freedom.
2. Define boundary conditions—applied displacements, applied tractions, concentrated
forces, body forces, initial strains.
3. Set up finite element mesh—define nodal coordinates, define element type, define ele-
ment connectivity (i.e. which nodes are attached to which elements).
4. Set up stiffness matrix—calculate local element stiffnesses and assemble global stiffness
matrix
5. Set up force matrix—calculate local element force and assemble global force matrix
6. Apply displacement boundary conditions—remove any known degrees of freedom from
the system of equations.
7. Invert stiffness matrix and solve for displacements.
8. Calculate strain field from displacements and stress field from strain field and consti-
tutive law.
9. Carry out postprocessing—plot deformed mesh, contours of stress strain etc.

Steps 1, 2, 9 and usually 3 are carried out by the user. Steps 4 to 8 and sometimes 3
(automatic meshing) are carried out by the finite element program.

181