Chapter 12

Three-dimensional Isoparametric Elements

Abstract Three-dimensional isoparametric elements are considered. Shape func-
tions for hexahedral elements with eight nodes and twenty nodes are given. A
method for computing derivatives of shape functions in the global coordinate system
is presented. Numerical integration is used for estimation of element matrices and
vectors. It is explained how to organize efcient computation of an element stiffness
matrix. Algorithms for calculating nodal equivalents of surface loads and stresses at
integration points and nodes are demonstrated.
12.1 Shape Functions
Hexahedral (or brick-type) linear eight-node and quadratic twenty-node three-
dimensional isoparametric elements are depicted in Figure 12.1. The term isopara-
metric means that the geometry and displacement elds are specied in parametric
form and are interpolated with the same functions. The shape functions used for in-
terpolation are polynomials of the local coordinates [, K and ] (1 [, K, ] 1).
Displacements are interpolated according to the following relations:
{u} = [N]{q},
{u} ={u v w},
{q} ={u
1
v
1
w
1
u
2
v
2
w
2
...},
(12.1)
where {u} is a displacement vector at a point with local coordinates [, K and ];
u, v, w are displacement components along global coordinate axes x, y and z; {q} is
a nodal displacement vector with entries u
i
, v
i
, w
i
(displacement values at nodes);
[N] is a matrix of shape functions.
Element shape is determined by coordinate interpolation using the matrix of
shape functions [N]:
129
1
3
4
2
1
2
3
4
5
6
7
8
5
6
7
8
9 10
11
12
13
14
15
16
17
18
19
20
(a) (b) (c)
x
h
z
Fig. 12.1 Linear (a) and quadratic (b) three-dimensional nite elements and their representation in
the local coordinate system (c)
{x} = [N]{x
e
},
{x} ={x y z},
{x
e
} ={x
1
y
1
z
1
x
2
y
2
z
2
...}.
(12.2)
Here, x, y, z are point coordinates and x
i
, y
i
, z
i
are coordinates of nodes. The matrix
of shape functions is dened as:
[N] =

N
1
0 0 N
2
0 0 ...
0 N
1
0 0 N
2
0 ...
0 0 N
1
0 0 N
2
...

. (12.3)
The shape functions of the three-dimensional linear element are:
N
i
=
1
8
(1 +[
0
)(1 +K
0
)(1 +]
0
),
[
0
= [[
i
,
K
0
= KK
i
,
]
0
=]]
i
.
(12.4)
A linear hexahedral element can be degenerated into a triangular prism (Fig-
ure 12.2a) by shrinking an element face into a straight line. Linear shape functions
are not affected by this degeneration.
For the quadratic element with twenty nodes, the shape functions can be written
in the following form:
1, 4
2
3
5, 8
6
7
1, 8, 7
2
3
4
6
5
9, 12
10
13, 20,19
11
14
15
17
18
16
(a) (b)
Fig. 12.2 Linear (a) and quadratic (b) elements degenerated into triangular prisms
N
i
=
1
8
(1 +[
0
)(1 +K
0
)(1 +]
0
)([
0
+K
0
+]
0
2) at vertices,
N
i
=
1
4
(1 [
2
)(1 +K
0
)(1 +]
0
) , i = 2, 6, 14, 18,
N
i
=
1
4
(1 K
2
)(1 +[
0
)(1 +]
0
) , i = 4, 8, 16, 20,
N
i
=
1
4
(1 ]
2
)(1 +[
0
)(1 +K
0
) , i = 9, 10, 11, 12.
(12.5)
In the above relations, [
i
, K
i
, ]
i
are values of local coordinates [, K, ] at nodes.
For the degenerated quadratic hexahedral element shown in Figure 12.2b six
shape functions for nodes located at a face opposite to the degenerated face should
be modied. In the case of the degeneration depicted in Figure 12.2b the shape
functions are modied in the following way [32]:
N

i
= N
i
+', i = 3, 5, 15, 17,
N

i
= N
i
2', i = 4, 16,
' =
1
16
_
1 [
2
__
1 K
2
_
(1 +]]
i
).
(12.6)
12.2 StrainDisplacement Matrix
The strain vector {H} contains six different components of the strain tensor:
{H} ={H
x
H
y
H
z
J
xy
J
yz
J
zx
}. (12.7)
The straindisplacement matrix for three-dimensional elements is:
[B] = [D][N] = [B
1
B
2
B
3
...], (12.8)
[B
i
] =

wN
i
wx
0 0
0
wN
i
wy
0
0 0
wN
i
wz
wN
i
wy
wN
i
wx
0
0
wN
i
wz
wN
i
wy
wN
i
wz
0
wN
i
wx

. (12.9)
The derivatives of three-dimensional shape functions with respect to global coordi-
nates are obtained as follows:

wN
i
wx
wN
i
wy
wN
i
wz

= [J]
1

wN
i
w[
wN
i
wK
wN
i
w]

, (12.10)
where the Jacobian matrix has the appearance
[J] =

wx
w[
wy
w[
wz
w[
wx
wK
wy
wK
wz
wK
wx
w]
wy
w]
wz
w]

. (12.11)
The partial derivatives of x, y, z with respect to [, K,] are found by differentiation
of displacements expressed through shape functions and nodal displacement values:
wx
w[
=

wN
i
w[
x
i
,
wx
wK
=

wN
i
wK
x
i
,
wx
w]
=

wN
i
w]
x
i
,
wy
w[
=

wN
i
w[
y
i
,
wy
wK
=

wN
i
wK
y
i
,
wy
w]
=

wN
i
w]
y
i
,
wz
w[
=

wN
i
w[
z
i
,
wz
wK
=

wN
i
wK
z
i
,
wz
w]
=

wN
i
w]
z
i
.
(12.12)
The transformation of integrals from the global coordinate system to the local coor-
dinate system is performed with the use of the determinant of the Jacobian matrix:
dV = dxdydz =|J|d[dKd]. (12.13)
12.3 Element Properties
Element matrices and vectors for the three-dimensional hexagonal element are given
by the following expressions:
stiffness matrix
[k] =
_
1
1
_
1
1
_
1
1
[B]
T
[E][B]|J|d[dKd], (12.14)
thermal vector (ctitious forces to simulate thermal expansion)
{h} =
_
1
1
_
1
1
_
1
1
[B]
T
[E]{H
t
}|J|d[dKd], (12.15)
force vector due to distributed face load
{p} =
_
1
1
_
1
1
[N]
T
{p
S
}
ds
d[
d[dK, (12.16)
equivalent stress vector (with negative sign)
{p
V
} =
_
1
1
_
1
1
_
1
1
[B]
T
{V}|J|d[dKd]. (12.17)
In the above relations [E] is the elasticity matrix of the material, {H
t
} is the vector
of thermal strains and {V} is the current stress vector. Three-time application of the
one-dimensional Gauss quadrature rule leads to the following numerical integration
procedure:
I =
_
1
1
_
1
1
_
1
1
f ([, K, ])d[dKd]
=
n

i=1
n

j=1
n

k=1
f ([
i
, K
j
, ]
k
)w
i
w
j
w
k
=
N

m=1
f ([
m
, K
m
, ]
m
)W
m
.
(12.18)
Here, [
i
, K
j
, ]
k
are the abscissas of numerical integrations; w
i
, w
j
, w
k
are the corre-
sponding weights. The third line of the above equation is a practical implementation
of three-dimensional integration rules as a single loop. In this case, N = n n n,
and W
m
are triple products of the Gauss weights.
Usually, 2 2 2 integration is used for three-dimensional linear elements and
333 integration is applied to evaluation of the stiffness matrix for quadratic el-
ements. For more efcient integration, a special 14-point Gauss-type rule [14, 17] is
employed, which provides sufcient precision of integration for three-dimensional
quadratic elements.
12.4 Efcient Evaluation of Element Matrices and Vectors
Calculation of the element stiffness matrix by multiplication of three matrices in-
volves many arithmetic operations with zeros. In addition, the three-dimensional
case does not contain any additional variants of the elasticity matrix and displace-
ment differentiation matrix as in the two-dimensional case. Because of this, a simple
and efcient three-dimensional algorithm for the stiffness matrix can be formulated
by performing matrix multiplications in closed form. Then, the coefcients of the
element stiffness matrix [k] are expressed as follows:
k
mn
ii
=
_
V
_
E
wN
m
wx
i
wN
n
wx
i
+P
_
wN
m
wx
i+1
wN
m
wx
i+1
+
wN
m
wx
i+2
wN
m
wx
i+2
__
dV,
k
mn
i j
=
_
V
_
O
wN
m
wx
i
wN
n
wx
j
+P
wN
m
wx
j
wN
n
wx
i
_
dV,
E = O +2P.
(12.19)
Here, O and P are elastic Lame constants; m, n are local node numbers; i, j are
indices related to coordinate axes (x
1
, x
2
, x
3
). The cyclic rule is employed in the
above equation if the coordinate indices exceed 3.
The analogous expression for the element thermal vector is
h
m
i
=
_
V
(3O +2P)DT
wN
m
wx
i
dV, (12.20)
where m is node number, i is coordinate index, D is the thermal-expansion coef-
cient, and T is temperature.
Aconvenient expression for estimation of the nodal vector equivalent to the stress
distribution is
(p
V
)
m
i
=
_
V
wN
m
wx
j
V
i j
dV. (12.21)
Here, we suppose summation over repeated index j. The nodal vector is calculated
with a negative sign for convenience of assembly.
12.5 Calculation of Nodal Equivalents for External Loads
The nodal equivalent of the external surface load is estimated according to Equation
3.20
{p} =
_
S
[N]
T
{p
S
}dS. (12.22)
g

x
h
x
h
Fig. 12.3 Determination of a normal vector g on the curved face of the nite element
Integration of the nodal equivalent of the surface load (12.16) requires explanation
related to expressing an area element dS through local coordinates [, K.
For integration of the surface load over the nite element face, it is necessary
to determine the outward normal vector at any point of the face. Vectors e
[
and
e
K
tangent to the local coordinates [ and K have the following components in the
global coordinate system (Figure 12.3):
e
[
=
_
wx
w[
wy
w[
wz
w[
_
,
e
K
=
_
wx
wK
wy
wK
wz
wK
_
.
(12.23)
Derivatives of the global coordinates with respect to local coordinates can be esti-
mated using two-dimensional shape functions:
wx
w[
=

wN
m
w[
x
m
, ...
wx
wK
=

wN
m
wK
x
m
, ... .
(12.24)
The normal vector g is equal to the vector product of the tangential vectors e
[
and
e
K
:
g = e
[
e
K
. (12.25)
The components of the outward normal vector g can be conveniently computed
through the following expressions:
g
x
=
wy
w[
wz
wK

wz
w[
wy
wK
,
g
y
=
wz
w[
wx
wK

wx
w[
wz
wK
,
g
z
=
wx
w[
wy
wK

wy
w[
wx
wK
.
(12.26)
Since integration of the surface load is performed in the local coordinate system [,
K, the Jacobian of transformation between the physical surface coordinates and the
local coordinates should be determined:
|J([, K)| =|g| =
_
g
2
x
+g
2
y
+g
2
z
. (12.27)
The element of area dS is expressed through local coordinates [, K as
dS =|J([, K)| d[dK, (12.28)
and the nodal equivalent for the surface load applied to a nite element face is
calculated in the form
[p] =
_
1
1
_
1
1
[N([, K)]
T
{p
S
([, K)} |J([, K)|d[dK. (12.29)
If the surface load is directed along the surface normal, the components of the unit
normal vector (normalized vector g) are used to produce components of the surface
load at integration points.
12.6 Example: Nodal Equivalents of a Distributed Load
Estimate the nodal equivalents of a distributed load with constant intensity p
s
= 1
applied to the at square face with size L = 1 of the three-dimensional quadratic
element shown in Figure 12.4.
Solution
The nodal equivalent of the distributed load can be calculated according to Equa-
tion 12.29. The calculation procedure is simplied if we take into account that the
element face is a at square with sides of length L = 1 parallel to coordinate axes x
and y. In this particular case
L
1 2 3
4
5
6 7
8
L
x
y
p
s
x
h
Fig. 12.4 Distributed load on a face of a three-dimensional quadratic element
dx
d[
=
1
2
,
dy
dK
=
1
2
and integration in (12.29) can be performed in a simpler way:
[p] =
_
1
1
_
1
1
[N]
T
{p
S
}
1
4
d[dK.
From symmetry properties, the equivalent nodal forces have identical values at all
corner nodes and at all midside nodes. It is sufcient to determine the nodal equiv-
alents at nodes 1 and 2.
Using two-dimensional shape functions [N
1
] and [N
2
] (10.7) nodal equivalents
for corner node 1 and midside node 2 are determined as:
p
1
=
_
1
1
_
1
1
1
4
(1 [)(1 K)(1 +[ +K)
1
4
d[dK =
1
12
,
p
2
=
_
1
1
_
1
1
1
2
(1 [
2
)(1 K)
1
4
d[dK =
1
3
.
Nodal equivalents of the distributed load with constant intensity for the face of the
three-dimensional quadratic element are shown in Figure 12.5.
At rst glance it seems unusual to observe nodal equivalents at corner nodes
directed in the opposite direction of the applied load. After counting the total force
12
1
12
1
12
1
12
1
3
1
3
1
3
1
3
1
Fig. 12.5 Distributed load on a face of a three-dimensional quadratic element
across the face, we nd that it is equal to unity conrming our calculations. Negative
directions of corner nodal forces are related to nonlinearity of shape functions.
12.7 Calculation of Strains and Stresses
After computing element matrices and vectors, the assembly process is used to com-
pose the global equation system. Solution of the global equation system provides
displacements at nodes of the nite element model. Using disassembly, nodal dis-
placement for each element can be obtained.
Strains inside an element are determined with the use of the displacement differ-
entiation matrix:
{H} = [B]{q}. (12.30)
Stresses are calculated with Hookes law:
{V} = [E]{H
e
} = [E]({H} {H
t
}), (12.31)
where {H
t
} is the vector of free thermal expansion:
{H
t
} ={DT DT DT 0 0 0}. (12.32)
It should be noted that displacement gradients (and hence strains and stresses) have
quite difference precision at different points inside nite elements. The highest pre-
cisions for displacement gradients are at the geometric center for the linear element
and at reduced integration points 22 2 for the quadratic hexagonal element.
12.8 Extrapolation of Strains and Stresses
For quadratic elements, displacement derivatives have the best precision at 222
integration points with local coordinates [, K, ] = 1/

3. In order to build a
continuous eld of strains or stresses, it is necessary to extrapolate result values from
2 2 2 integration points to the vertices of a twenty-node element (numbering of
integration points and vertices is shown in Figure 12.6).
Results are calculated at 8 integration points, and a trilinear extrapolation in the
local coordinate system [, K, ] is used:
f
i
= L
i(m)
f
(m)
, (12.33)
where f
(m)
are known function values at reduced integration points, f
i
are function
values at vertex nodes, and L
i(m)
is the extrapolation matrix:
1
2
3 4
( ) 3
(1)
( ) 2
( ) 4
5
6
7
8
( ) 5
( ) 6
( ) 7
( ) 8
Fig. 12.6 Numbering of integration points and vertices for the twenty-node isoparametric element
L
i(m)
=

A B B C B C C D
B C C D A B B C
C D B C B C A B
B C A B C D B C
B A C B C B D C
C B D C B A C B
D C C B C B B A
C B B A D C C B

,
A =
5 +

3
4
, B =

3+1
4
,
C =

31
4
, D =
5

3
4
.
(12.34)
Stresses are extrapolated fromintegration points to all nodes of elements. Values for
midside nodes can be calculated as an average between values for two vertex nodal
values. Then averaging of contributions from the neighboring nite elements is per-
formed for all nodes of the nite element model. Averaging produces a continuous
eld of secondary results specied at nodes of the model with quadratic variation
inside nite elements. Later, the results can be interpolated to any point inside an
element or on its surface using quadratic shape functions.
Problems
12.1. Derive explicit expressions for shape functions N
6
and N
15
for the three-
dimensional quadratic element. Use the general relations (12.5).
12.2. Show that three-dimensional linear shape functions (12.4) satisfy the equality