Week 10 MECH3361

1
Recap Finite Element Method (FEM)
Bar Element Method 1 (Displacement Method) Method 2 (Energy Method)
Displacement u
j i
u L x u L x x u ) / ( ) / 1 ( ) ( + =
j i j j i i
u u u N u N u + = + = ) 1 ( ) ( ) ( ) (
c
L L
u u
i j
A
=

= c
( ) Bu u N Nu =
|
.
|

\
|
= = =
dx
d
dx
d
dx
du
c
o
L
E
E
A
= = c o Bu E E = = c o
Force or W=U A = A
|
.
|

\
|
=
A
= = k
L
EA
A
L
E
A F o ( ) f u u B B u
T
2
1
2
1
=
(
(

}
V
T T
dV E U
Equilibrium eqn f ku= f ku=
Stiffness matrix
(

=
(

=
1 1
1 1
L
EA
k k
k k
k
(

=
1 1
1 1
L
EA
k

Finite Element (FE) equilibrium equation: f ku=
1D
Spring
element

)
`

=
)
`

j
i
j
i
f
f
u
u
k k
k k


1D Bar
Elements
:

)
`

=
)
`

j
i
j
i
f
f
u
u
L
EA
1 1
1 1



2D Bar
Element

(
(
(
(

0
0
'
'
'
'
'
'
2 2
2 2
2 2
2 2
j
i
j
j
i
i
f
f
v
u
v
u
m lm m lm
lm l lm l
m lm m lm
lm l lm l
L
EA



Beam
Element



Distributed Load Equivalent nodal load
Single
bar
element

Multiple
bar
elements



Beam
element



Week 10 MECH3361
2
8.5 Two-Dimensional Problems

General Formula for the Stiffness Matrix
Recall the displacement in terms of shape function in bar element
Nu =
)
`

=
j
i
j i
u
u
N N u ] [
We can extend it to 2D, where the displacements (u, v) in a plane element are interpolated
from nodal displacements (u
i
, v
i
) using shape functions N
i
as follows,

+ + = =
+ + = =

2 2 1 1
2 2 1 1
) , ( ) , ( ) , (
) , ( ) , ( ) , (
v y x N v y x N y x v v
u y x N u y x N y x u u

In the matrix form: Nd u =

=
)
`

2
2
1
1
2 1
2 1
0 0
0 0
v
u
v
u
N N
N N
v
u

where N is the shape function matrix, u the displacement vector inside each element and d
the elemental nodal displacement vector. Here we assume that u depends on the nodal values
of u only, and v on nodal values of v only.
From strain-displacement relation, the strain vector can be derived as:
( )d DN Du = = c
Or Bd = c
where ( ) DN B= is the strain-displacement matrix and D is derivative matrix (see below).

Consider the strain energy stored in an element and take the 2D as an example:
} } }
} }
= = ) ( =
+ + = =
V
T
V
T
V
T
V
xy xy yy yy xx xx
V
T
dV dV dV
dV dV U
) )
2
1
2
1
2
1
) (
2
1
] [
2
1
d d E( ( Ec c c Ec
c o c o c o c o

Since nodal displacement vector d is independent on the elemental coordinate dV.
kd d d d
T
V
T T
dV U =
|
|
.
|

\
|
=
}
E
2
1

From this, we obtain the general formula for the element stiffness matrix,

Remarks:
- Note that unlike 1-D cases, E here is a matrix that is given by the stress-strain relation
(generalised Hookes rule)
- The stiffness matrix k is symmetric since E is symmetric.
- Also note that for given the material property E, the behaviour of k depends on the B
matrix only, which in turn on the shape functions. Thus, the quality of finite elements
in representing the behaviour of a structure is entirely determined by the choice of
shape functions.
Week 10 MECH3361
3
Constant Strain Triangle (CST or T3)
This is the simplest 2-D element, which is also called linear triangular element.
For this element, we have three nodes at the vertices of
the triangle, which are numbered around the element in
the counter-clockwise direction. Each node has two
degrees of freedom (can move in both x and y
directions). The displacements u and v are assumed to be
linear functions within the element, that is,

y b x b b y x v v
y b x b b y x u u
6 5 4
3 2 1
) , (
) , (
+ + = =
+ + = =

where b
i
(i = 1, 2, ..., 6) are constants to be determined.
From these, the strains are found to be,

( )
2 3 2 1
b y b x b b
x x
u
xx
= + +
c
c
=
c
c
= c


( )
6 6 5 4
b y b x b b
y y
v
yy
= + +
c
c
=
c
c
= c

( ) ( )
5 3 3 2 1 6 5 4
2 b b y b x b b
y
y b x b b
x x
v
y
u
xy
+ = + +
c
c
+ + +
c
c
=
|
|
.
|

\
|
c
c
+
c
c
= c
Obvious the strains are constant throughout the element. Thus, we have the name constant
strain triangle (CST) element.
Displacement functions should satisfy the following six equations,

3 6 3 5 4 3 3 3
3 3 3 2 1 3 3 3
2 6 2 5 4 2 2 2
2 3 2 2 1 2 2 2
1 6 1 5 4 1 1 1
1 3 1 2 1 1 1 1
) , (
) , (
) , (
) , (
) , (
) , (
y b x b b y y x x v v
y b x b b y y x x u u
y b x b b y y x x v v
y b x b b y y x x u u
y b x b b y y x x v v
y b x b b y y x x u u
+ + = = = =
+ + = = = =
+ + = = = =
+ + = = = =
+ + = = = =
+ + = = = =

Solving these six equations, we can find the coefficients b
1
, b
2
, ..., b
6
in terms of nodal
displacements and coordinates. Substituting these coefficients into displacement equations
and rearranging the terms, we obtain:

=
)
`

3
3
2
2
1
1
3 2 1
3 2 1
0 0 0
0 0 0
v
u
v
u
v
u
N N N
N N N
v
u
(linear distribution)
where the shape functions (note that they are linear functions of x and y) are
| |
| |
| | y x x x y y y x y x
A
N
y x x x y y y x y x
A
N
y x x x y y y x y x
A
N
) ( ) ( ) (
2
1
) ( ) ( ) (
2
1
) ( ) ( ) (
2
1
1 2 2 1 1 2 2 1 3
3 1 1 3 3 1 1 3 2
2 3 3 2 2 3 3 2 1
+ + =
+ + =
+ + =

Week 10 MECH3361
4
where
(
(
(

=
3 3
2 2
1 1
1
1
1
det
2
1
y x
y x
y x
A is the area of the triangle
Using the strain-displacement relation, we have,

(
(
(

= =

3
3
2
2
1
1
12 21 31 13 23 32
21 13 32
12 31 23
0 0 0
0 0 0
2
1
v
u
v
u
v
u
y x y x y x
x x x
y y y
A
xy
yy
xx
Bd
c
c
c

where x
ij
= x
i
- x
j
and y
ij
= y
i
- y
j
(i, j = 1, 2, 3). Again, we see constant strains within the
element (because variables x and y disappeared). From stress-strain relation (Hookes law),
we see that stresses obtained in the CST element should also be constant.

Applying general formula for stiffness matrix, we obtain it for the CST element,

in which t is the thickness of the element. Notice that k for CST is a 6 by 6 symmetric matrix.
The matrix multiplication can be carried out by a computer program.

Applications of the CST Element:
- Use in the areas where the strain gradient is small.
- Use in the mesh transition areas (fine mesh to coarse mesh) to accommodate different
mesh topologies.
- Avoid using CST in stress concentration or other crucial areas in the structure, such as
edges of holes and corners.
- Recommended for quick and preliminary FE analysis of 2-D problems.

Linear Strain Triangle (LST or T6)
This element is also called quadratic triangular element.

There are six nodes on this element: three corner nodes
and three mid-side nodes. Each node has two degrees of
freedom (DOF) as before. The displacements (u, v) are
assumed to be quadratic functions of (x, y),

2
12 11
2
10 9 8 7
2
6 5
2
4 3 2 1
y b xy b x b y b x b b v
y b xy b x b y b x b b u
+ + + + + =
+ + + + + =

where b
i
(i=1,2,12) are the constants to be determined.

Thus the strains are found to be:
( ) y b x b b y b xy b x b y b x b b
x x
u
xx 5 4 2
2
6 5
2
4 3 2 1
2 + + = + + + + +
c
c
=
c
c
= c
( ) y b x b b y b xy b x b y b x b b
y y
v
yy 12 11 9
2
12 11
2
10 9 8 7
2 + + = + + + + +
c
c
=
c
c
= c
Week 10 MECH3361
5

( )
( )
( ) ( ) ( )y b b x b b b b
y b xy b x b y b x b b
y
y b xy b x b y b x b b
x x
v
y
u
xy
11 6 10 3 5 3
2
12 11
2
10 9 8 7
2
6 5
2
4 3 2 1
2 2 2
2
+ + + + + =
+ + + + +
c
c
+
+ + + + +
c
c
=
|
|
.
|

\
|
c
c
+
c
c
= c

Obviously all these strains all are in linear functions. Thus, we have the linear strain
triangle (LST), which may provide better accuracy than the CST.

In a natural coordinate system, the six shape functions for the LST element can be defined:
| |
q q q q
q q q q
) 1 ( 4 ), 1 ( 4 , 4
1 ) 1 ( 2 ) 1 ( ), 1 2 ( ), 1 2 (
6 5 4
3 2 1
= = =
= = =
N N N
N N N

Displacement can be written as:
i
i
i i
i
i
v N v u N u

= =
= =
6
1
6
1
, (quadratic distribution)
which means that for given nodal displacements, one can know the displacement at any point
within the element. In other words, we use shape functions and nodal displacement to express
displacement field (function).

Remarks: Note that the element stiffness matrix is still given by dV
V
T e
B E B k
}
= , but
EB B
T
is quadratic w.r.t. x and y. In general, the integral has to be computed numerically.

Bi-Linear Quadrilateral Element (Q4)

There are four nodes at the corners of the
quadrilateral shape. In the natural coordinate system
(, q), the four shape functions are given as,
) 1 )( 1 (
4
1
1
q = N , ) 1 )( 1 (
4
1
2
q + = N ,
) 1 )( 1 (
4
1
3
q + + = N , ) 1 )( 1 (
4
1
4
q + = N
Mapping from Cartesian coordinate to natural
coordinate systems can be done as follows:

i
i
i i
i
i
y N y x N x

= =
= =
4
1
4
1
) , ( , ) , ( q q ,
Displacement can be written as:
i
i
i i
i
i
v N v u N u ) , ( , ) , (
4
1
4
1
q q

= =
= = (bilinear distribution)
Strain:
|
|
.
|

\
|
c
c
=
c
c
=

=
i
i
i xx
u N
x x
u
) , (
6
1
q c

|
|
.
|

\
|
c
c
=
c
c
=

=
i
i
i yy
v N
y y
v
6
1
) , ( q c

|
|
.
|

\
|
c
c
+
|
|
.
|

\
|
c
c
=
|
|
.
|

\
|
c
c
+
c
c
=

= =
i
i
i i
i
i xy
u N
y
v N
x x
v
y
u
6
1
6
1
) , ( ) , ( 2 q q c
We cannot directly obtain the results because u and v are defined as functions of natural
coordinate and q rather than as a function of x and y. We can however use chain rule:
Week 10 MECH3361
6

c
c
c
c
(
(
(
(

c
c
c
c
c
c
c
c
=

c
c
c
c
y
u
x
u
y x
y x
u
u
J

] [
q q

q

or

c
c
c
c
(

c
c
c
c
(
(
(
(

c
c
c
c
c
c
c
c
=

c
c
c
c

q q

u
u
J J
J J
u
u
y x
y x
y
u
x
u
J
*
22
*
21
*
12
*
11
1
] [
1


where [J] is called J acobian matrix.

Coefficients in [J] can be obtained as:


|
|
.
|

\
|
c
c
=
c
c
|
|
.
|

\
|
c
c
=
c
c
i
i
i
i
y
N y
x
N x

,
Thus
q
c
c
c
+
c
c
=
c
c
=
u
J
u
J
x
u
x
*
12
*
11

where
*
12
*
11
, J J are the coefficients in the first row of the inverse Jacobian matrix [J]
-1
.


|
|
.
|

\
|
c
c
=
c
c
|
|
.
|

\
|
c
c
=
c
c
i
i
i
i
u
N u
u
N u
q q
,

Thus we change the integration from Cartesian coordinate to natural coordinate as

} } }

= =
1
1
1
1
q d d t dV
T
V
T e
J EB B B E B k
which needs to use numerical integration method.

According to Gaussian rule of integration, one can have:
( ) ( )

} }
= =

= =
n
i
m
j
j i j i
W W d d I
1 1
1
1
1
1
, , q | q q |
Gaussian integration points:

Gaussian points and weights:
m,n q, W
i

11 0.0
22 5774 . 0 3 / 3 =
1
33
0
7746 . 0 6 . 0 =
0.8888
0.5555
The stress, strain results at the Gaussian point should be most accuratein the entire element.
3
3
=
3
3
= q
3
3
=
3
3
= q

q
x
y
Week 10 MECH3361
7
Quadratic Quadrilateral Element (Q8)

) 1 )( 1 )( 1 (
4
1
) 1 )( 1 )( 1 (
4
1
) 1 )( 1 )( 1 (
4
1
) 1 )( 1 )( 1 (
4
1
4
3
2
1
q q
q q
q q
q q
+ + =
+ + =
+ + =
+ + =
N
N
N
N

) 1 )( 1 (
2
1
) 1 )( 1 (
2
1
) 1 )( 1 (
2
1
) 1 )( 1 (
2
1
2
8
2
7
2
6
2
5
q
q
q
q
=
+ =
+ =
=
N
N
N
N

Mapping from Cartesian coordinate to natural coordinate systems:
i
i
i i
i
i
y N y x N x

= =
= =
8
1
8
1
, ,
The displacement field is given by:
i
i
i
i
i
v N v u N u

= =
= =
8
1
8
1
, (quadratic distributions)
The stiffness matrix:

} } }

= =
1
1
1
1
q d d t dV
T
V
T e
J EB B B E B k
which needs to use numerical integration method as
( ) ( )

} }
= =

= =
n
i
m
j
j i j i
W W d d I
1 1
1
1
1
1
, , q | q q |

Remarks:
- Q4 and T3 are usually used together in a mesh with linear elements.
- Q8 and T6 are usually applied in a mesh composed of quadratic elements.
- Quadratic elements are preferred for stress analysis, because of their high accuracy
and the flexibility in modeling complex geometry, such as curved boundaries.

Example 8.7 Sketch the displacement distribution of the following elements (CST, LST, Q4
and Q8) and explain the difference (final exam 2010, 2011, 2012):

Solution:
CST element (T3): linear distribution of displacement and constant in the strain.

4-nodel element (Q4): bilinear distribution of displacements and bilateral linear in the strain

q
x
y
x
y

q
mapping
1
2
3
4
5
6
7
8
1
2
3
5
6
7
8
4
Week 10 MECH3361
8

3-node CST element 4-nodel bilinear element

LST element (6-node): quadratic distribution of displacement and linear distribution in the strain

8-node element: quadratic distribution of displacements and linear distribution in the strain

6-node LST element 8-nodel quadratic element


Example 8.8
For the 8-node element in Example 8.7, if 22 Gaussian points are considered, go on to
determine the displacements of all the Gaussian points if the nodal displacement vectors are
given as below (exam 2010, 2012):

(

=
)
`

2 1 0 1 1 1 1 2
1 2 1 0 1 1 0 1
v
u


Soln:
The displacement functions as per shape functions N
i
and nodal displacements u
i
and v
i
can be
expressed as:

i
i
i i
i
i
v N v u N u ) , ( , ) , (
4
1
4
1
q q

= =
= =
Shape functions for 8-node elements can be written in terms of natural coordinate:

) 1 )( 1 )( 1 (
4
1
) , (
) 1 )( 1 )( 1 (
4
1
) , (
) 1 )( 1 )( 1 (
4
1
) , (
) 1 )( 1 )( 1 (
4
1
) , (
4
3
2
1
q q q
q q q
q q q
q q q
+ + =
+ + =
+ + =
+ + =
N
N
N
N

) 1 )( 1 (
2
1
) , (
) 1 )( 1 (
2
1
) , (
) 1 )( 1 (
2
1
) , (
) 1 )( 1 (
2
1
) , (
2
8
2
7
2
6
2
5
q q
q q
q q
q q
=
+ =
+ =
=
N
N
N
N

x
y
u
i
j
k
u
j
u
i
u
k x
y
i
j
k
v
j
v
i
v
k
v
v
l
l
x
y
u
i
j
k
u
j
u
i
u
k x
y
i
j
k
v
j
v
i
v
k
v
v
l
l
Quadratic surface
Quadratic surface
Week 10 MECH3361
9
Gaussian points can be defined in the natural coordinates as (refer to Table on Page 6)
Gaussian point #1: =0.5774, q=0.5774; Gaussian point #2: =-0.5774, q=0.5774;
Gaussian point #3: =-0.5774, q=-0.5774; Gaussian point #4: =0.5774, q=-0.5774;

Therefore for { }
T
u 1 2 1 0 1 1 0 1 } { = :
6923 . 1 ) 5774 . 0 , 5774 . 0 (
) 1 ( ) 5774 . 0 1 )( 5774 . 0 1 ( 5 . 0 ) 2 ( ) 5774 . 0 1 )( 5774 . 0 1 ( 5 . 0
) 1 ( ) 5774 . 0 1 )( 5774 . 0 1 ( 5 . 0 ) 0 ( ) 5774 . 0 1 )( 5774 . 0 1 ( 5 . 0
) 1 ( ) 5774 . 0 5774 . 0 1 )( 5774 . 0 1 )( 5774 . 0 1 ( 25 . 0
) 1 ( ) 5774 . 0 5774 . 0 1 )( 5774 . 0 1 )( 5774 . 0 1 ( 25 . 0
0 ) 5774 . 0 5774 . 0 1 )( 5774 . 0 1 )( 5774 . 0 1 ( 25 . 0
1 ) 5774 . 0 5774 . 0 1 )( 5774 . 0 1 )( 5774 . 0 1 ( 25 . 0
) 1 )( 1 ( 5 . 0
) 1 )( 1 ( 5 . 0 ) 1 )( 1 ( 5 . 0 ) 1 )( 1 ( 5 . 0 ) 1 )( 1 )( 1 ( 25 . 0
) 1 )( 1 )( 1 ( 25 . 0 ) 1 )( 1 )( 1 ( 25 . 0 ) 1 )( 1 )( 1 ( 25 . 0
) 5774 . 0 , 5774 . 0 ( ) , (
2 2
2 2
8
2
7
2
6
2
5
2
4
3 2 1
8 8 7 7 6 6 5 5 4 4 3 3 2 2 1 1
8
1
=
+ + +
+ + +
+ +
+ +
+ +
+ + =
+
+ + + + + + +
+ + + + + + =
+ + + + + + + = = =

=
u
u
u u u u
u u u
u N u N u N u N u N u N u N u N u N u u
i
i
i
q
q q q q q
q q q q q q
q

For { }
T
v 2 1 0 1 1 1 1 2 } { = :
9553 . 0 ) 5774 . 0 , 5774 . 0 (
) 2 ( ) 5774 . 0 1 )( 5774 . 0 1 ( 5 . 0 ) 1 ( ) 5774 . 0 1 )( 5774 . 0 1 ( 5 . 0
) 0 ( ) 5774 . 0 1 )( 5774 . 0 1 ( 5 . 0 ) 1 ( ) 5774 . 0 1 )( 5774 . 0 1 ( 5 . 0
) 1 ( ) 5774 . 0 5774 . 0 1 )( 5774 . 0 1 )( 5774 . 0 1 ( 25 . 0
) 1 ( ) 5774 . 0 5774 . 0 1 )( 5774 . 0 1 )( 5774 . 0 1 ( 25 . 0
1 ) 5774 . 0 5774 . 0 1 )( 5774 . 0 1 )( 5774 . 0 1 ( 25 . 0
2 ) 5774 . 0 5774 . 0 1 )( 5774 . 0 1 )( 5774 . 0 1 ( 25 . 0
) 5774 . 0 , 5774 . 0 ( ) , (
2 2
2 2
8 8 7 7 6 6 5 5 4 4 3 3 2 2 1 1
8
1
=
+ + +
+ + +
+ +
+ +
+ +
+ + =
+ + + + + + + = = =

=
v
v N v N v N v N v N v N v N v N v N v v
i
i
i
q

Similarly, displacements in other Gaussian points can be found using Matlab code as:
xi=0.5774; eta=0.5774;
u1=1; u2=0; u3=-1; u4=-1; u5=0; u6=1; u7=2; u8=1;
v1=2; v2=1; v3=-1; v4=-1; v5=1; v6=0; v7=-1; v8=-2;
N1=-(1-xi)*(1-eta)*(1+xi+eta)/4; N2=-(1+xi)*(1-eta)*(1-xi+eta)/4;
N3=-(1+xi)*(1+eta)*(1-xi-eta)/4; N4=-(1-xi)*(1+eta)*(1+xi-eta)/4;
N5=(1-xi*xi)*(1-eta)/2; N6=(1+xi)*(1-eta*eta)/2;
N7=(1-xi*xi)*(1+eta)/2; N8=(1-xi)*(1-eta*eta)/2;
u=N1*u1+N2*u2+N3*u3+N4*u4+N5*u5+N6*u6+N7*u7+N8*u8
v=N1*v1+N2*v2+N3*v3+N4*v4+N5*v5+N6*v6+N7*v7+N8*v8
Use the Matlab code, we can find the following results:
6923 . 1 ) 5774 . 0 , 5774 . 0 ( = u , 9553 . 0 ) 5774 . 0 , 5774 . 0 ( = v
6218 . 1 ) 5774 . 0 , 5774 . 0 ( = u , 7956 . 1 ) 5774 . 0 , 5774 . 0 ( = v
3075 . 1 ) 5774 . 0 , 5774 . 0 ( = u , 3778 . 0 ) 5774 . 0 , 5774 . 0 ( = v
0445 . 1 ) 5774 . 0 , 5774 . 0 ( = u , 1290 . 0 ) 5774 . 0 , 5774 . 0 ( = v