Transcript The Finite Element Method A Practical Course
The Finite Element Method
A Practical Course
CHAPTER 5: FEM FOR 2D SOLIDS
CONTENTS
INTRODUCTION LINEAR TRIANGULAR ELEMENTS – Field variable interpolation – Shape functions construction – – Using area coordinates Strain matrix – Element matrices LINEAR RECTANGULAR ELEMENTS – – – – – Shape functions construction Strain matrix Element matrices Gauss integration Evaluation of
m
e
CONTENTS
LINEAR QUADRILATERAL ELEMENTS – – – – Coordinate mapping Strain matrix Element matrices Remarks HIGHER ORDER ELEMENTS COMMENTS (GAUSS INTEGRATION) CASE STUDY
INTRODUCTION
2D solid elements are applicable for the analysis of plane strain and plane stress problems.
A 2D solid element can have a rectangular or quadrilateral shape with curved edges.
triangular , straight or 2D solid element can deform only in the plane of the 2D solid.
At any point, there are two components in
x
and directions for the displacement as well as forces.
y
INTRODUCTION
For plane strain problems, the thickness of the element is unit , but for plane stress problems, the actual thickness must be used.
In this course, it is assumed that the element has a uniform thickness
h
.
Formulating 2-D elements with given variation of thickness is also straightforward, as the procedure is the same as that for a uniform element.
2D solids – plane stress and plane strain Plane stress Plane strain
LINEAR TRIANGULAR ELEMENTS
Less accurate than quadrilateral elements Used by most mesh generators for complex geometry Linear triangular element Triangular elements Nodes
y, v
1 (
x
1
, y
1 ) (
u
1 ,
v
1 ) 3 (
x
3
, y
3 ) (
u
3 ,
v
3 )
f sy A f sx
2 (
x
2
, y
2 ) (
u
2 ,
v
2 )
x, u
Field variable interpolation
U
h
N
( , )
d
e
where
N
N
1 0
d
e
u v u v
u v
2 3 1 1 2 3 displaceme displaceme displaceme nts nts at node 2 nts at node at node 1 3
y, v
0
N
1
N
0 2 0
N
2
N
0 3 0
N
3 Node 1 Node 2 Node 3 (Shape functions) 1 (
x
1
, y
1 ) (
u
1 ,
v
1 ) 3 (
x
3
, y
3 ) (
u
3 ,
v
3 )
f sy A f sx
2 (
x
2
, y
2 ) (
u
2 ,
v
2 )
x, u
Shape functions construction
Assume,
N
1 1 1 1
N
2
a
2 2 2
N
3 3 3 3
N i
i i
i
or
i
= 1, 2, 3
N i
1
x
p
T y
a
b i
p
T
Shape functions construction
Delta function property:
i
( ,
j
Solving,
j
) Therefore, 1 for 0 for
i i
j j
1 ( , 1 1 1 ( , 2 2 ) 0 3 ) 0 1 1 ( , 1 1 ) 1 ( , 2 2 ) 1 1
b x
1 1
b x
1 2
c y
1 1
c y
1 2 1 0 ( , 3 3 ) 1
b x
1 3
c y
1 3 0
a
1
x y
2 3
x y
3 2 ,
b
1 2
A e
y
2
y
3 ,
c
1 2
A e x
3
x
2 2
A e
Shape functions construction
A e
1 2
P
1 1 1 2 1
x
1
x
2
x
3
y
1
y
2
y
3 1 2 [(
x y
2 3
x y
3 2 Area of triangle Moment matrix
y
2
y x
3 ) 1 (
x
3 2 ) ] 1 Substitute
a
1 ,
b
1 and
c
1 back into
N
1
= a
1
+ b
1
x + c
1
y
:
N
1 1 2
A e
[(
y
2
y
3 )(
x
x
2
x
3
x
2 )(
y
y
2 )]
Shape functions construction
Similarly, 2 ( , 1 1 ) 0 2 ( , 2 2 2 ( , 3 3 ) 0 3 ( , 2 1 ) 0 2 ) 0 3 ( , 3 3
N
2 1 2
A e
[(
x y
3 1
x y
1 3 1 2
A e
[(
y
3
y
1 )(
x
x
3
y
3
y x
1 ) (
x
1
x y
3 ) ]
x
1
x
3 )(
y
y
3 )]
N
3 1 2
A e
[(
x y
1 2
x y
1 1 1 2
A e
[(
y
1
y
2 )(
x
x
1
y
1
y x
2 ) (
x
2
x y
1 ) ]
x
2
x
1 )(
y
y
1 )]
Shape functions construction
N
i
i i
i
where
a i
1 2
A e
(
x y j k
x y k j
)
b i
1 2
A e
(
y j
y k
)
c i
1 2
A e
(
x k
x j
)
i J
= 1, 2, 3 ,
k
permutation
k
= 3, 1
k
determined from cyclic
i i
= 1, 2
j j
= 2, 3
Using area coordinates
Alternative method of constructing shape functions
y i,
1 2-3-P:
P k
, 3
A
1
A
1 1 1 1 2 1
x x
2
y y
2 1 2 [(
x y
2 3
x y
3 2
x
3
y
3
L
1
A
1
A e
Similarly, 3-1-P
A
2
y
2
y x
3 ) (
x
3
x y
2 ) ]
L
2
A
2 1-2-P
A
3
A e A
3
j
, 2
L
3
A e x
Using area coordinates
Partitions of unity:
L
1
L
2
L
3 1
L
1
L
2
L
3
A
1
A e A
2
A e
A
3
A e
A
1
A
2
A
3
A e
1 Delta function property: e.g.
L
1 = 0 at if
P
at nodes 2 or 3 Therefore,
N
1
L
1 ,
N
2
L
2 ,
N
3
L
3
U
h
N
( , )
d
e
Strain matrix
xx
yy
xy
x v
u y
u y
v x
LU
LU B
LN
x
0
y
LNd
e
Bd
e
0
y
x
N
where
L
x
0
y
0
y
x
B
b
1 0
c
1 0
c
1
b
1
b
2 0
c
2 0
c
2
b
2
b
3 0
c
3 0
c
3
b
3 (constant strain element)
Element matrices
k
e
V e T
B cB
d
V
A e
0
h T
B cB
d
A
A e h T
B cB
d
A
Constant matrix
k
e
hA e T
B cB m
e
V e
T
N N
d
V
A e
0
h
d
z
T
N N
d
A
A e h
T
N N
d
A
Element matrices
For elements with uniform density and thickness,
m
e
h
A e
N N
1 0
N N
2 1 0
N N
3 1 0 1 0
N N
1 1 0 0
N N
3 1 1
N N
1 2 0
N N
2 2 0
N N
3 2 0 0
N N
1 2 0 0
N N
3 2 2
N N
1 3 0
N N
2 3 0
N N
3 3 0 0
N N
1 3 0 0
N N
3 3 3 d
A
Eisenberg and Malvern (1973):
A L m
1
L n
2
L
3
p
d
A
(
m
m
!
n
!
p
!
n
p
2 )!
2
A
Element matrices
m
e
hA
12 2 0 2
sy
.
1 0 2 0 1 0 2 1 0 1 0 2 0 1 0 1 0 2
y, v
1 (
x
1
, y
1 ) (
u
1 ,
v
1 ) 3 (
x
3
, y
3 ) (
u
3 ,
v
3 )
f sy A f sx
2 (
x
2
, y
2 ) (
u
2 ,
v
2 )
x, u
f
e
l
N
T
f f
d
l
Uniform distributed load:
f
e
1 2
l
2 3
f f f
0 0
f y x y x
LINEAR RECTANGULAR ELEMENTS
Non-constant strain matrix More accurate representation of stress and strain Regular shape makes formulation easy
Shape functions construction
Consider a rectangular element
d
e u
u
1 2
v u v u
2 3 3 displacements at node 1 displacements at node 2 displacements at node 3 displacements at node 4
y, v
4 (
x
4 ,
y
4 ) (
u
4 ,
v
4 )
2b
1 (
x
1 ,
y
1 ) (
u
1 ,
v
1 )
2a
3 (
x
3 ,
y
3 ) (
u
3 ,
v
3 )
f sy f sx
2 (
x
2 ,
y
2 ) (
u
2 ,
v
2 )
x, u
Shape functions construction
4 ( 1, +1) (
u
4 ,
v
4 )
y, v
4 (
x
4 ,
y
4 ) (
u
4 ,
v
4 )
2b
1 (
x
1 ,
y
1 ) (
u
1 ,
v
1 )
2a
3 (
x
3 ,
y
3 ) (
u
3 ,
v
3 )
f sy f sx
2 (
x
2 ,
y
2 ) (
u
2 ,
v
2 )
2
1 ( 1, 1)
2
(
u
1 ,
v
1 )
x, u
x
(
x
2
a
x
1 ) / 2 ,
y
(
y
2
b
y
1 ) / 2 3 (1, +1) (
u
3 ,
v
3 ) 2 (1, 1) (
u
2 ,
v
2 )
U
h
N
( , )
d
e
where
N
N
1 0 Node 1 0
N
1
N
2 0 Node 2 0
N
2
N
3 0 Node 3 0
N
3
N
4 0 0
N
4 Node 4
Shape functions construction
N
1
N
2
N
3
N
4 1 4 ( 1 )( 1 ) 1 4 ( 1 )( 1 ) 1 4 ( 1 )( 1 ) 1 4 ( 1 )( 1 )
N j
1 4 ( 1
j
)( 1
j
) 4 ( 1, +1) (
u
4 ,
v
4 )
N
3 at node 1
N
3 at node 2
N
3 at node 3
N
3 at node 4 1 4 (1 )(1 ) 1 1 1 4 (1 )(1 ) 1 1 0 0 1 4 (1 )(1 ) 1 1 1 1 4 (1 )(1 ) 1 1 0 Delta function property
2
1 ( 1, 1)
2
(
u
1 ,
v
1 )
i
4 1
N i
N
1 1 4 [(1 )(1 1 4 [2(1
N
2
N
3
N
4 )(1 )(1 )(1 )] Partition of unity 3 (1, +1) (
u
3 ,
v
3 ) 2 (1, 1) (
u
2 ,
v
2 )
Strain matrix
B
LN
4 1
a
0 1
b
0 1 1
b
a
1
a
0 1
b
0 1 1
b
a
1
a
0 1
b
0 1 1
b
a
1
a
0 1
b
0 1 1
b
a
Note: No longer a constant matrix!
Element matrices
x a
,
y b
d
x
d
y
=
ab
d d Therefore,
k
e
A
h
B
T
cB
d
A
1 1 1 1
ab h
B
T
cB
d d
m
e
V
T
N N
d
V
A
0
h
d
z
T
N N
d
A
A
h
T
N N
d
A
1 1 1 1
abh
T
N N
Element matrices
f
e
l
N
T
f f
d
l
For uniformly distributed load,
f
e
b
0 0
f f f f
0 0
x y x y
y, v 2b
4 (
x
4 ,
y
4 ) (
u
4 ,
v
4 ) 1 (
x
1 ,
y
1 ) (
u
1 ,
v
1 )
2a
3 (
x
3 ,
y
3 ) (
u
3 ,
v
3 )
f sy f sx
2 (
x
2 ,
y
2 ) (
u
2 ,
v
2 )
x, u
Gauss integration
For evaluation of integrals in
k
e
and
m
e
(in practice) In 1 direction:
I
1 1
f
( )d
j m
1
w j f
(
j
)
m
gauss points gives exact solution of polynomial integrand of
n
= 2
m
- 1 In 2 directions:
I
1 1 1 1
f
i n x n
1
j y
1
w w f i j
i j
)
Gauss integration
m Gauss Point
j Gauss Weight w j
3 4 1 2 5 6 0 -1/ 3, 1/ 3 0.6, 0, 0.6
2 1, 1 5/9, 8/9, 5/9 -0.861136, -0.339981, 0.339981, 0.861136 -0.906180, -0.538469, 0, 0.538469, 0.906180 -0.932470, -0.661209, -0.238619, 0.238619, 0.661209, 0.932470 0.347855, 0.652145, 0.652145, 0.347855 0.236927, 0.478629, 0.568889, 0.478629, 0.236927 0.171324, 0.360762, 0.467914, 0.467914, 0.360762, 0.171324
Accuracy order
n
1 3 5 7 9 11
Evaluation of
m
e
m
e
hab
9 4 0 4
sy
.
2 0 4 0 2 0 4 1 0 2 0 4 2 0 4 0 1 0 0 2 0 4 2 0 1 0 2 0 4 0 2 0 1
Evaluation of
m
e m ij
hab
1 1 1 1
hab
16
hab
1 1 ( 1 ( 1 1
N i
4 3
i
N
i
j j d
d
)( 1 )( 1 1 3
j
i
j
)
d
) 1 1 ( 1
i
)( 1
j
)
d
E.g.
m
33
hab
( 1 4 1 3 1 1 )( 1 1 3 1 1 ) 4
hab
9 Note: In practice, gauss integration is often used
LINEAR QUADRILATERAL ELEMENTS
Rectangular elements have limited application Quadrilateral elements with unparallel edges are more useful Irregular shape requires coordinate mapping before using Gauss integration
Coordinate mapping
y
4 (
x
4 ,
y
4 ) 3 (
x
3 ,
y
3 ) 4 ( 1, +1) 3 (1, +1) 1 (
x
1 ,
y
1 ) 2 (
x
2 ,
y
2 )
x
1 ( 1, 1) 2 (1, 1) Physical coordinates Natural coordinates
U
h
N
d
e
X
N
x
e
(Interpolation of displacements) (Interpolation of coordinates)
Coordinate mapping
X
N
x
e
where
X
N
1
N
2
N
3
N
4 1 4 ( 1 )( 1 ) 1 4 ( 1 )( 1 ) 1 4 ( 1 )( 1 ) 1 4 ( 1 )( 1 ) ,
x
e x y x
1 1
y
y x y
2 3 3 4 4 coordinate at node 1 coordinate at node 2 coordinate at node 3 coordinate at node 4
x
i
4 1
N i
( , )
x i y
i
4 1
N i
( , )
y i
Coordinate mapping
Substitute 1 into
x
i
4 1
x y
1 2 1 2 (1 )
x
2 (1 )
y
2 1 2 1 2 (1 )
x
3 (1 )
y
3 or
N i
( , )
x i x
y
1 2 (
x
2 1 2 (
y
2
x
3 )
y
3 ) 1 2 (
x
3 1 2 (
y
3
x
2
y
2 ) ) Eliminating ,
y
(
y
3 (
x
3
y
2
x
2 ) ) {
x
1 2 (
x
2
x
3 )} 1 2 (
y
2
y
3 )
y
4 (
x
4 ,
y
4 ) 3 (
x
3 ,
y
3 ) 4 ( 1, +1) 3 (1, +1) 1 (
x
1 ,
y
1 ) 2 (
x
2 ,
y
2 )
x
1 ( 1, 1) 2 (1, 1)
Strain matrix
N i
N
i
N
x i
N
x i
x
x
N i
y
N i
y
y
y
or
N i
N
i
J
N i
x
N i
y
where
J
x
x
Since
X
y
y
(Jacobian matrix)
N
x
e
,
J
N
1
N
1
N
2
N
2
N
3
N
3
N
4
N
4
x x
2
x
1 3
x
4
y y y
3 1 2
y
4
Strain matrix
Therefore,
N i
x
N i
y
J
1
N
i
N
i
(Relationship between differentials of shape functions w.r.t. physical coordinates and differentials w.r.t. natural coordinates) Replace differentials of
N i
with differentials of
N i
w.r.t.
x
w.r.t. and and
y
B
LN
0
x
y
0
y
x
N
Element matrices
Murnaghan (1951) : d
A
=det |
J
| d d
k
e
1 1 1 1
h
T
B cB
det
J m
e
V
N
T
N
d
V
A
0
h
d
x
N
T
N
d
A
1 1 1 1
h
N
T
N
det
J
d d
A
h
N
T
N
d
A
Remarks
Shape functions used for interpolating the coordinates are the same as the shape functions used for the interpolation of the displacement field. Therefore, the element is called
isoparametric element
.
Note that the shape functions for coordinate interpolation and displacement interpolation do not have to be the same.
Using the different shape functions for coordinate interpolation and displacement interpolation, respectively, will lead to the development of so-called
superparametric
elements.
subparametric
or
HIGHER ORDER ELEMENTS
(
p
,0,0) Higher order triangular elements (0,0,
p
)
n d
= (
p
+1)(
p
+2)/2 (1,0,
p
1) (0,1,
p
1) Node
i
,
I
(2,0,
p
2) Argyris, 1968 :
L
1
L
2 (
p
1,1,0)
i
(
I
,
J
,
K
)
L
3
K p N i
I I l L l
1
J J L l
2 )
K K
(
L
3 ) (0,
p
1,1)
l
(
L
) ( (
L
L
L
0 )(
L
L
0 )(
L
L
1 )
L
1 ) (
L
(
L
L
1) )
L
( 1) ) (0,
p
,0)
HIGHER ORDER ELEMENTS
Higher order triangular elements (Cont’d)
y, v
3
N
1
N
4 (2
L
1 4 2 1)
L
1 5 6 2 1 4 Quadratic element
x, u y, v
1 9 8 3 7 4 10 5 6
N
1 1 2 (3
L
1 1)(3
L
1 2)
L
1
N
4
N
10 9 2
L L
1 2 27 (3
L
1
L L L
1 2 3 1) 2
x, u
Cubic element
HIGHER ORDER ELEMENTS
(0,
m
) Higher order rectangular elements Lagrange type: (Zienkiewicz et al., 2000) (
n
,
m
)
N i
1
D N N
1
D J
l I n
l J m
i
(
I
,
J
) 0
l k n
( 0 )( 1 (
k
0 )(
k
1 ) ) ( (
k k k
1 1 )( )(
k
k
1
k
) 1 ) ( (
k
n
)
n
) (0,0) (
n
,0)
HIGHER ORDER ELEMENTS
Higher order rectangular elements(Cont’d) 4 7 3 J=2 J=1 8 9 6 (9 node quadratic element) J=0 1 I=0
N
1
N
0 1
D N
2
N
2 1
D N
3
N
2 1
D N
4
N
0 1
D N
0 1
D N
0 1
D N
2 1
D N
2 1
D
5 I=1 2 I=2 1 4 (1 1 4 (1 1 4 (1 )(1 1 4 (1 )(1 ) )
N
5
N
1 1
D N
6
N
2 1
D N
7
N
1 1
D N
8
N
0 1
D N
9
N
1 1
D N
0 1
D N
1
N
1 2
D
1
D N
1 1
D N
1 1
D
1 2 (1 )(1 )(1 1 2 (1 )(1 )(1 ) 1 2 (1 )(1 )(1 1 2 (1 )(1 2 )(1 2 )
HIGHER ORDER ELEMENTS
8 Higher order rectangular elements(Cont’d) Serendipity type: 4 =1 7 3 0 6
N j N j N j
1 4 (1
j
)(1
j j
j
1)
j
1, 2, 3, 4 1 2 (1 2 )(1
j
)
j
5, 7 1 2 (1
j
)(1 2 )
j
6, 8 1 = 1 5 2 (eight node quadratic element)
HIGHER ORDER ELEMENTS
Higher order rectangular elements(Cont’d) 4 10 9 1 11 12 5 6 2 3 8 7
N j
1 32 (1
j
)(1
j
2 9 2 10) for corner nodes
j
1, 2, 3, 4
N j
9 32 (1
j
)(1 2 for side nodes
j
j
) 7, 8, 11, 12 where
j N j
9 32 (1
j
)(1 2 )(1 9
j
) for side nodes
j
5, 6, 9, 10 where
j
1 and 1 3 and
j j
1 3 1 (twelve node cubic element)
ELEMENT WITH CURVED EDGES
3 3 5 5 6 6 2 2 4 4 1 1 4 1 8 7 5 3 2 6 8 1 4 7 3 6 5 2
COMMENTS (GAUSS INTEGRATION)
When the Gauss integration scheme is used, one has to decide how many Gauss points should be used.
Theoretically, for a one-dimensional integral, using
m
points can give the exact solution for the integral of a polynomial integrand of up to an order of (2
m
1).
As a general rule of thumb, more points should be used for higher order of elements.
COMMENTS (GAUSS INTEGRATION)
Using smaller number of Gauss points tends to counteract the
over-stiff behaviour
associated with the displacement based method.
Displacement in an element is assumed using shape functions. This implies that the deformation of the element is somehow prescribed in a fashion of the shape function.
This prescription gives a constraint to the element. The so constrained element behaves stiffer than it should be. It is often observed that higher order elements are usually softer than lower order ones. This is because using higher order elements gives less constraint to the elements.
COMMENTS (GAUSS INTEGRATION)
Two gauss points for linear elements, and two or three points for quadratic elements in each direction should be sufficient for many cases. Most of the explicit FEM codes based on explicit formulation tend to use one-point integration to achieve the best performance in saving CPU time.
CASE STUDY
Side drive micro-motor
10N/m
CASE STUDY
10N/m 10N/m Elastic Properties of Polysilicon Young’s Modulus, E Poisson’s ratio, Density, 169GPa 0.262
2300kgm -3
CASE STUDY
Analysis no. 1: Von Mises stress distribution using 24 bilinear quadrilateral elements (41 nodes)
CASE STUDY
Analysis no. 2: Von Mises stress distribution using 96 bilinear quadrilateral elements (129 nodes)
CASE STUDY
Analysis no. 3: Von Mises stress distribution using 144 bilinear quadrilateral elements (185 nodes)
CASE STUDY
Analysis no. 4: Von Mises stress distribution using 24 eight-nodal, quadratic elements (105 nodes)
CASE STUDY
Analysis no. 5: Von Mises stress distribution using 192 three-nodal, triangular elements (129 nodes)
Analysis no.
1 2 3 4 5
CASE STUDY
Number / type of elements Total number of nodes in model Maximum Von Mises Stress (GPa) 24 bilinear, quadrilateral 41 0.0139
96 bilinear, quadrilateral 129 0.0180
144 bilinear, quadrilateral 24 quadratic, quadrilateral 192 linear, triangular 185 105 129 0.0197
0.0191
0.0167