18 April 2017 1
49325 Computer-Aided Mechanical Design
Lecture 04: FEM for Two-Dimensional Problems
Isoparametric Elements and their Numerical
Integration
Dr. Zhen (Jeff) Luo
School of EMMS, FEIT, UTS
Office: CB11.09.132
Phone: 9514 2994;
Email: [email protected]
The term “Isoparametric”, which means “same parameter”, is
derived from the use of the same shape functions [N] to define
the shape/coordinates as are used to define the displacements of
any point within the element.
x
y
x y 1 1 ,
x y 2 2 , x y 3 3 ,
x y 4 4 ,
1
2
3
4
u v 1 1 ,
u v 2 2 ,
u v 3 3 ,
u v 3 3 , { } and ( , ) ? X X N(ξ,η)
{ } and ( , ) ? d u N(ξ, η)
1: “Isoparametric”?How to map the local coordinates for any point in the square
reference element to the global coordinates for the quadrilateral
physical element ? For instance, considering the square reference
element, the node 1 has , while in the quadrilateral physical
element we have .
x y
1, 1
x x y y 1 1 ,
x
y
x y 1 1 ,
x y 2 2 ,
x y 3 3 ,
x y 4 4 ,
1
2
3
4
1, 1 1, 1
1,1 1,1
1 2
4 3
Square reference element ? Quadrilateral physical element
Reference: FEM Meshing:
Local Coordinate System Global Coordinate System4
(1) For isoparametric finite elements: The
displacement field, the strain and the stress
are all to be formulated in
Local Coordinate System
(2) However, please note that strain, stress
and strain are, originally, all defined in
Global Coordinate System
(3) Hence, a transformation has to be applied
to mapping or project quantities between
Local - Global Coordinate Systems
1 2
3 4
1 1
1 1 , 1 1
4 4
1 1
1 1 , 1 1
4 4
N N
N N
1 1 2
1 2 3 4 2
1 2 3 4 3
Shape Functions: ( , )
3 4 4
FEM Meshing
( , ) ( , ) 0 0 0 0
( , ) ( , ) 0 0 0 0
i i
i i
x y x
x N x y N N N N
y N y N N N N x
y x y
N
Shape functions: Element shapes
Global Coordinate System
Local Coordinate System
Eq. (1):
Local Coordinate System
x x ( , ) ?{ }
x
y
x y 1 1 , x y 2 2 ,
x y 3 3 ,
1
2
3
4
( , ) x y
x y 4 4 ,
o
o
Step 1: Select type/shape of
isoparametric element6
“Isoparametric” – Two identical shape functions to define both the
coordinates and the displacements of any point within an element.
The same shape functions, used to compute the coordinates of any point
(x, y) within the element, is also used to calculate the displacements of this
point within the element:
Step 2: Select displacement functions
1 1 2
1 2 3 4 2
1 2 3 4 3
3 4 4
( , ) ( , ) 0 0 0 0
( , )
( , ) ( , ) 0 0 0 0
i i
i i
u v u
u N u v N N N N
v N v u N N N N
v u v
N(ξ,η)
FEM:{d}
u
u d ( , ) ?{ }
1 1 2
1 2 3
FEM Mes
4 2
1 2 3 4 3
4 4
n
3
hi g
( , ) ( , ) 0 0 0 0
,
( , ) ( , ) 0 0 0 0
i i
i i
x y x
x N x y N N N N
y N y N N N N x
y x y
N(ξ,η)
X
Ref to:
Eq. (2):
Eq. (1):
( , ) ( , ) { }
FEM Mesh
x x N
( , ) ( , ) { }
FEM Mesh
u d N7
Step 3: Define Strain-Displacement Relationship
For 2D problems, the three strains (two normal strains; one engineering shear
strain) corresponds to the three stresses. The strains are originally defined in
terms of the derivatives of the displacements u and v with respect to global x
and y coordinates. However, (x, y) are also the functions of . We will need
to express the derivatives of the displacement in the global (x, y) coordinate
system in terms of its derivatives in local coordinate system:
1 0 0 0
( ) 0 0 0 1
0 1 1 0
,
x y
xy
uxuy
uxvy
u v
y x v
vx y
T
{ε }
,
Eq. (3) ?
( , ) ( , ) ( , ) ?{ } u v u v u v x y x y
FEM
ε ε ε d
,
,
,
?
,8
22 12
21 11
22 12
21 11
0 0
1 0 0
det 0 0
0 0
uxu vxv
u
J J u
J J
J J v
J J
y y
v
Α
J
22 12
21 11
1
det
u
J J
J J
uxuy
u
J
22 12
21 11
1
det
v
J J
J J
vxvy
v
J
Then the first-order derivatives for the displacements of any point (u,v)
in the element with respect to the Global Coordinates (x,y) are given as
! ?
1
1 2 4 3
1 2
1 2 4 3
2 3
1 2 4 3
3
1 2 4 3
0 0 0 0
0 0 0 0
0 0 0 0
0 0 0 0
u N N N N u
v u
u N N N N
v u
v N N N N
v
v N N N N u
S
4 4
v
{d}
1 1 2
1 2 3 4 2
1 2 3 4 3
3 4 4
, ,
,
,
( ) ( ) 0 0 0 0
( )
( , ) ( ) 0 0 0 0
i i
i i
u v u
u N u v N N N N
v N v u N N N N
v u v
N(ξ,η)
FEM:{d}
u
Recall
Eq. (2):
!Known Information given by Shape Functi
(1 ) 0 (1 ) 0 (1 ) 0 (1 ) 0
1 (1 ) 0 (1 ) 0 (1 ) 0 (1 ) 0
4 0 (1 ) 0 (1 ) 0 (1 ) 0 (1 )
0 (1 ) 0 (1 ) 0 (1 ) 0 (1 )
u u v v
1 1 2 2
ons:
3 3 4 4
u v u v u v u v
S
FEM:{d}
31 2 4 3
1 2 4
1
(1 ) (1 ) (1 ) (1 )
4 1
(1 ) (1 ) (1 ) (1 )
4
NN N N N
N N N
1 2
3 4
1 1
1 1 , 1 1
4 4
1 1
1 1 , 1 1
4 4
N N
N N
Substitute into [S]:
!
Global Coordinate System
Local Coordinate System
Local Coordinate Systemwhere the matrix [B] is called the strain-displacement matrix, it is only
dependent on nodal local coordinates of the element, as well as the
derivatives of shape functions. The more detailed information about [B]:
, ,
B
ε T A ξ, η S ξ, η d = B d
In this case, in the numerical implementation of FEM, we can calculate
the strain-displacement matrix in the local/natural system:
B B B B B , , , , , 3 8 1 2 3 4 3 2 3 2 3 2 3 2
22 12
3 2 21 11
21 11 22 12
0
1
( ) 0 : 1, 2, 4
e
, 3,
d t
i i
i i
i
i i i i
Where
N N
J J
N N
J J nodes i
N N N N
J J J J
B
J18/04/2017 12
1 1
22 12
1 1
1 21 11
1 1 1 1
21 11 22 12
( ) ( )
0
1 ( ) ( )
( ) 0
det
( ) ( ) ( ) ( )
, ,
, ,
,
, , , ,
N N
J J
N N
J J
N N N N
J J J J
B
J
2 2
22 12
2 2
2 21 11
2 2 2 2
21 11 22 12
( ) ( )
0
1 ( ) ( )
( ) 0
det
( ) ( ) ( ) ( )
, ,
, ,
,
, , , ,
N N
J J
N N
J J
N N N N
J J J J
B
J
3 3
22 12
3 3
3 21 11
3 3 3 3
21 11 22 12
( ) ( )
0
1 ( ) ( )
( ) 0
det
( ) ( ) ( ) ( )
, ,
, ,
,
, , , ,
N N
J J
N N
J J
N N N N
J J J J
B
J
4 4
22 12
4 4
4 21 11
4 4 4 4
21 11 22 12
( ) ( )
0
1 ( ) ( )
( ) 0
det
( ) ( ) ( ) ( )
, ,
, ,
,
, , , ,
N N
J J
N N
J J
N N N N
J J J J
B
JIn this case, in the numerical implementation of FEM, we can calculate
the strain within the local (or natural) coordinate system:
( , ) ?{ }
FEM
ε d
1 0 0 0
( ) 0 0 0 1
0 1 1 0
,
x y
xy
uxuy
uxvy
u v
y x v
vx y
T
{ε }
Eq. (3):
22 12
21 11
22 12
21 11
0 0
1 0 0
det 0 0
0 0
uxu vxv
u
J J u
J J
J J v
J J
y y
v
Α
J
Known Information given by Shape Functi
(1 ) 0 (1 ) 0 (1 ) 0 (1 ) 0
1 (1 ) 0 (1 ) 0 (1 ) 0 (1 ) 0
4 0 (1 ) 0 (1 ) 0 (1 ) 0 (1 )
0 (1 ) 0 (1 ) 0 (1 ) 0 (1 )
u u v v
1 1 2 2
ons:
3 3 4 4
u v u v u v u v
S
FEM:{d}
!
!
,
,
!
, ,
B
ε T A ξ,η S ξ,η d = B d
For linear elastic materials, the stress-strain relations come from the generalized
Hooke’s law. Particularly, for isotropic materials, there are only two independent
material properties: Young’s Modulus E and Poisson’s Ratio v. Consider a
differential element in the body, the Hook’s law gives :
14
2 2
(1 )
(1 )
2(1 )
x x y
y x y
xy xy xy
E
v
v
E
v
v
E
G
v
,
2
1 0
1 0
1
0 0 1 2
x x
y y
xy xy
v
E
v
v
v
σ D ε
(4) Stress-Strain Relationship
,
, , ,
ε
σ D ε B D d
Hence, the stress of isotropic materials is
( , ) ? { }
FEM
σ d
Eq. (4):15
Unlike the Constant-Strain-Triangular
(CST) element, the strain and stress in
the quadrilateral Q4 element are not
constant within the element. Instead,
they are functions of and
consequently vary within the element.
( , ) 16
(5) Derive Element Stiffness Matrix and EquationsKeep in mind that for the basic plane stress element, the total potential
energy is now a function of the nodal displacement {d}, such that the
partial derivative of ∏ with respect to the nodal displacement {d} is
T 0
V
dV
{d} [B] D [B] {d} f
0
e
e
T
V
e
dV
k
[B] D [B] {d} f
So, for element with a constant thickness h in the global x-y system, we have
17
( , ) ( , )
e
e T
A
k [B ] D [B ] x y x y hdxdy
e
e T
V
k [B] D [B] dV
Note: The strain energy is originally defined in
Global Coordinate System18
1, 1 1, 1
1,1 1,1
1 2
4 3
k [B ] D [B ] J e T 8 8 1 1 1 1 ( , ) ( , ) ( , ) 8 3 3 8 3 3 h d d
How to calculate [k] in ξ-η
( , ) ( , )
e
e T
A
k [B ] D [B ] x y x y hdxdy
dxdy d d J
Local Coordinate System:
Global Coordinate System:18 April 2017 19
Multidimensional Gauss rules are obtained by successive application
of 1D Gauss rules. In two dimensions,
1 1 1 1 1 1 , , , i i i j i j
i j i
I y d d W y d WW y
In two dimensions, for example, a four-node Gauss rule (a 2-by-2 rule:
displacement field in x or y direction, linear distribution, 2 points will be good
enough for fully representing the linear distribution) is shown in the
following Figures, with i=1,2 and j=1,2 yields
I y d d W W y W W y W W y W W y 1 1 1 1 , , , , , 1 1 1 1 2 1 2 1 1 2 1 2 2 2 2 2
1 1
2 1
1 2
2 2
Gauss point 1: , 0.5773, 0.5773
Gauss point 2: , 0.5773, 0.5773
Gauss point 3: , 0.5773,0.5773
Gauss point 4: , 0.5773,0.5773
1 2
1.000
Weights:
1.000
W W
3. Numerical integration
(1) Gauss quadrature for stiffness matrix:20
1
2
3 4
1 1 , 2 1 ,
1 2 , 2 2 ,
j 2
j 1
i 1 i 2
0.5773 0.5773
0.5773
0.5773
Figure: Gaussian quadrature in two dimension (4 computational points)Numerical integration for K e
For an isoparametric plane element, the element stiffness can be computed by
using the numerical integration (e.g. Gaussian Quadrature), stated by
k [B ] D [B ] J e T 8 8 1 1 1 1 ( , ) ( , ) ( , ) 8 3 3 8 3 3 h d d
Local Coordinate System:
(1) For i=1 to 4
End
1 1 1 1 1 1 1 1
2 2 2 2 2 2 2 2
3 3 3 3 3 3 3 3
4 4 4 4 4 4 4 4
( , ) ( , ) ( , )
( , ) ( , ) ( , )
( , ) ( , ) ( , )
( , ) ( , ) ( , )
e T
T T T
hWW
hW W
hW W
hW W
k [B ] D [B ] J
[B ] D [B ] J
[B ] D [B ] J
[B ] D [B ] J
4
1
e T ( , ) ( , ) ( , ) i i i i i i i i
i
hWW
k [B ] D [B ] J (2) For i=1 to 2
For j=1 to 2
8 8
1 1 8 3 1 1 3 8 1 1 1 1
2 1 8 3 2 1 3 8 2 1 2 1
1 2 8 3 1 2 3 8 1
( , ) ( , ) ( , )
( , ) ( , ) ( , )
With respect to the Element Stiffness Matrix:
( , ) ( , ) ( ,
e
T T T
hWW
hW W
k
[B ] D [B ] J
[B ] D [B ] J
[B ] D [B ] J
2 1 2
2 2 8 3 2 2 3 8 2 2 2 2
)
( , ) ( , ) ( , ) T
hWW
[B ] D [B ] J hW W
2 2
1 1
e T ( , ) ( , ) ( , ) i j i j i j i j
j i
hWW
k [B ] D [B ] J
End (loop j=1,2)
End (loop i=1,2)23
Numerical Case:
(how to calculate [K] in local coordinate?)
1 1
2 1
1 2
2 2
Gauss point 1: , 0.5773, 0.5773
Gauss point 2: , 0.5773, 0.5773
Gauss point 3: , 0.5773,0.5773
Gauss point 4: , 0.5773,0.5773
1 2
1.000
1.000
W W
1 2
3 4
1 1
1 1 , 1 1 ,
4 4
1 1
1 1 , 1
, ,
, , 1
4 4
i j i j i j i
i
j
j i j i j i j
N N
N N
31 2 4 3
1 2 4
1
(1 ) (1 ) (1 ) (1 ) ;
4 1
(1 ) (1 ) (1 ) (1 ) ;
,
4
,i j i
j j j j
i i
j
i i
N N N N N
N N N N N
% Loop the 4 Gauss integration points (ξi, ηj)
For i=1 to 2
For j=1 to 2
(1) Calculate Shape functions with respect to local coordinate (ξ, η)
(2) Calculate first derivatives of shape functions with respect to (ξ, η)
1 1
11 12 2 2
21 22 3 3
4 4
1 1 1 1 1
4 1 1
,
1 1
j j j j
i j
i i i i
x y
J J x y
J J x y
x y
J
(3): Calculation of the Jacobian matrix J in (ξ, η)(4): Calculation the determinant of the Jacobian matrix |J| in (ξ, η)
1 2
1 2 3 4
3 4
0 1 1
1 1 0 1
det
8 1 0 1
1 1 0
, ( , )
j j i i
j i
i
i j
i j i j
i i j j
j i j
y y
x x x x
y y
J J
B B B B B i j i j i j i j i j , , , , , 1 2 3 4
(5): Calculation [B] in (ξ, η) based on (2), (3) and (4)
1 1
22 12
1 1
1 21 11
1 1 1 1
21 11 22 12
( , ,
, ,
) ( )
0
1 ( ) ( )
( ) 0
det
( ) ( ) ( ) ( )
,
,
, , , ,
i j i j
i j i j
i j
i j
i j i j i j i j
N N
J J
N N
J J
N N N N
J J J J
B
J
2 2
22 12
2 2
2 21 11
2 2 2 2
21 11 22 12
( , ,
, ,
) ( )
0
1 ( ) ( )
( ) 0
det
( ) ( ) ( ) ( )
,
,
, , , ,
i j i j
i j i j
i j
i j
i j i j i j i j
N N
J J
N N
J J
N N N N
J J J J
B
J26
End (loop j=1,2)
End (loop i=1,2)
(6): Calculation [K] based on (5)
1 1 1 1 1 1 1 1
2 1 2 1 2 1 2 1
1 2 1 2 1 2 1 2
2 2 2 2 2 2 2 2
( , ) ( , ) ( , )
( , ) ( , ) ( , )
( , ) ( , ) ( , )
( , ) ( , ) ( , )
e T
T T T
hWW
hW W
hWW
hW W
k [B ] D [B ] J
[B ] D [B ] J
[B ] D [B ] J
[B ] D [B ] J
k [B ] D [B ] J ( , ) ( , ) ( , ) i j i j i j i j T hW W
3 3
22 12
3 3
3 21 11
3 3 3 3
21 11 22 12
( ) , ,
, ,
,
,
,
( )
,
0
1 ( ) ( )
( ) 0
det
( ) ( ) ( ) ( ) , ,
i j i j
i j i j
i j
i j
i j i j i j
N N
J J
N N
J J
N N N N
J J J J
B
J
4 4
22 12
4 4
4 21 11
4 4 4 4
21 11 22 12
( , ,
, ,
) ( )
0
1 ( ) ( )
( ) 0
det
( ) ( ) ( ) ( )
,
,
, , , ,
i j i j
i j i j
i j
i j
i j i j i j i j
N N
J J
N N
J J
N N N N
J J J J
B
J27
In Q4 element, the stresses are not constant within the quadrilateral element,
because [B] is a function of , and the stress is also a function of . In
general, the stresses are evaluated at the same Gauss points as used to
evaluate the stiffness matrix (for 2-by-2 integration, 4 sets of stress data). To
reduce the data, it is often practical to evaluate the stresses at .
Element Q4 has the defect of displaying spurious shear stress when bending.
In practice, shear stress in Q4 element should be calculated at the center of
the element, the element stiffness can be computed at the center ,
which happens to be the Gauss point location of one-point quadrature rule.
It is well known that the Q4 element typically generates overstiffness finite
element results. So, for the Q4 element, an obvious benefit of one-point
integration is that the resulting [K] is not stiffened by spurious shear strain
because bending modes produce zero shear strain at the element center.
, ,
0
0
( , )
( 0, 0)
( 0, 0)
( 0, 0)
σ = D B d
BJ J
Element Strain and Stress:
(2) Gauss quadrature for stress calculation:28
1
2
3 4
1 1 , 2 1 ,
1 2 , 2 2 ,
j 2
j 1
i 1 i 2
Figure: Gaussian quadrature in two dimension
(central computational point)
0, 0