Transcript MANE 4240 & CIVL 4240 Introduction to Finite Elements Prof. Suvranu De Mapped element geometries and shape functions: the isoparametric formulation.

MANE 4240 & CIVL 4240
Introduction to Finite Elements
Prof. Suvranu De
Mapped element geometries
and shape functions: the
isoparametric formulation
Reading assignment:
Chapter 10.1-10.3, 10.6 + Lecture notes
Summary:
• Concept of isoparametric mapping
• 1D isoparametric mapping
• Element matrices and vectors in 1D
• 2D isoparametric mapping : rectangular parent elements
• 2D isoparametric mapping : triangular parent elements
• Element matrices and vectors in 2D
For complex geometries
© 2002 Brooks/Cole Publishing / Thomson Learning™
General quadrilateral elements
Elements with curved sides
Consider a special 4-noded rectangle in its local coordinate
system (s,t)
t
2
1
1
1
Displacement interpolation
u  N1u1  N 2 u 2  N 3u3  N 4 u 4
1
s
1
3
4
v  N1v1  N 2 v2  N 3 v3  N 4 v4
Shape functions in local coord system
1
N 1 ( s, t )  (1  s )(1  t )
4
1
N 2 ( s, t )  (1  s )(1  t )
4
1
N 3 ( s, t )  (1  s )(1  t )
4
1
N 4 ( s, t )  (1  s )(1  t )
4
Recall that
N1  N 2  N 3  N 4  1
Rigid body modes
N1 s1  N 2 s 2  N 3 s3  N 4 s 4  s
N1t1  N 2 t 2  N 3t 3  N 4 t 4  t
Constant
strain states
Goal is to map this element from local coords to a general
t
quadrilateral element in global coord system
2
t
2
1
1
x  x ( s, t )
y
1
1
1
y  y ( s, t )
3
s
s
1
s  s ( x, y )
s
4
4
3
t  t ( x, y )
x
Local coordinate system
Global coordinate system
N1 (s, t )
N1 ( x, y)
N1 (s, t )  N1 (s( x, y),t ( x, y))  N1 ( x, y)
In the mapped coordinates, the shape functions need to satisfy
1. Kronecker delta property
 1 at nodei
Ni  
0 at all other nodes
Then
2. Polynomial completeness
N
i
1
i
N x
i
i
x
i
N
i
i
yi  y
The relationship
x   N i ( s, t ) x i
i
y   N i ( s, t ) y i
i
Provides the required mapping from the local coordinate system
To the global coordinate system and is known as isoparametric
mapping
(s,t): isoparametric coordinates
(x,y): global coordinates
Examples
t
t
1
1
1
s
1
y
x
t
s
t
1
s
1
s
y
x
1D isoparametric mapping
3 noded (quadratic) element
1
3
1
2
1
Local (isoparametric)
coordinates
s 1  s 
N1 (s)  
2
s 1  s 
N 2 (s) 
2
N 3 ( s)  1  s 2
s
x1
x3 x2
1
3
x
2
Isoparametric mapping
3
x   N i ( s ) xi
i 1
s1  s 
s1  s 
x
x1 
x 2  1  s 2 x3
2
2


NOTES
1. Given a point in the isoparametric coordinates, I can obtain the
corresponding mapped point in the global coordinates using the
isoparametric mapping equation
s 1  s 
s 1  s 
x
x1 
x 2  1  s 2 x3
2
2

At s  1; x  x1
At s  0; x  x3
At s  1; x  x2
Question
x=? at s=0.5?

2. The shape functions themselves get mapped
In the isoparametric coordinates (s) they are polynomials.
In the global coordinates (x) they are in general nonpolynomials
Lets consider the following numerical example
x1  0; x2  6; x3  4
4
1
Isoparametric mapping x(s)
x

s 1  s 
2
s 1  s 
x1 
0
2
 4  3s  s 2
s 1  s 
2
s 1  s 
2

2
3
x
2

x2  1  s 2 x3


6  1  s2 4
Simple polynomial
Inverse mapping s(x)
s
3  25  4 x
2
Complicated function
Now lets compute the shape functions in the global coordinates
s1  s 
2
1  3  25  4 x  3  25  4 x 
1 

 



2
2
2


1
 10  x  2 25  4 x
2
 N 2 ( x)
N 2 (s) 


Now lets compute the shape functions in the global coordinates
N2(x)
N2(s)
1
3
1
2
1
1
4
s
N2(s) is a simple polynomial
s 1  s 
N 2 (s) 
2
1
1
2
3
x
2
N2(x) is a complicated function
N 2 ( x) 

1
10  x  2 25  4 x
2
However, thanks to isoparametric mapping, we always ensure
1. Knonecker delta property
2. Rigid body and constant strain states

Element matrices and vectors for a mapped 1D bar element
1
3
1
2
1
x
s
1
3
2
Displacement interpolation u  N1u1  N 2u2  N3u3  N d
Strain-displacement relation  
Stress   E  E Bd
dN
dN
du dN1

u1  2 u 2  3 u 3  B d
dx dx
dx
dx
 dN1 dN2 dN3 
The strain-displacement matrix B  

 dx dx dx 
The only difference from before is that the shape functions are in the
isoparametric coordinates
s 1  s 
2
s 1  s 
N 2 (s) 
2
N 3 ( s)  1  s 2
N1 (s)  
We know the isoparametric mapping
3
x   N i ( s ) xi
i 1
And we will not try to obtain explicitly the inverse map.
How to compute the B matrix?
Using chain rule
dNi ( s) dNi ( s) ds

dx
ds dx
Do I know
dN i ( s )
?
ds
Do I know
ds
?
dx
(*)
3
x   N i ( s ) xi
I know
i 1
3
Hence dx   dNi (s) xi  J ( Jacobianof m apping)
ds
From (*)
i 1
ds
dNi ( s)
dx
1 dNi ( s)

J ds
What does the Jacobian do?
dx
 J ds
Maps a differential element from the isoparametric coordinates to
the global coordinates
The strain-displacement matrix
 dN1 dN 2 dN3 
B

 dx dx dx 
1  dN1 dN 2 dN3 
 
J  ds ds ds 
For the 3-noded element
dNi ( s)
2s  1
2s  1
J 
xi 
x1 
x2  2sx3
ds
2
2
i 1
3
1  2s  1 2s  1

B 
 2s 
J 2
2

The element stiffness matrix
x2
k   EA B B dx
T
x1
1
  EA B B Jds  dx  Jds
T
1
NOTES
1. The integral on ANY element in the global coordinates in now an
integral from -1 to 1 in the local coodinates
2. The jacobian is a function of ‘s’ in general and enters the integral.
The specific form of ‘J’ is determined by the values of x1, x2 and x3.
Gaussian quadrature is used to evaluate the stiffness matrix
3. In general B is a vector of rational functions in ‘s’
Isoparametric mapping in 2D: Rectangular parent elements
© 2002 Brooks/Cole Publishing / Thomson Learning™
Parent element
Isoparametric mapping
Mapped element in
global coordinates
x   N i ( s, t ) x i
i
y   N i ( s, t ) y i
i
Shape functions of parent element in isoparametric
coordinates
1
(1  s )(1  t )
4
1
N 2 ( s, t )  (1  s )(1  t )
4
1
N 3 ( s, t )  (1  s )(1  t )
4
1
N 4 ( s, t )  (1  s )(1  t )
4
N1 ( s, t ) 
Isoparametric mapping
x   N i ( s, t ) x i
i
y   N i ( s, t ) y i
i
NOTES:
1. The isoparametric mapping provides the map (s,t) to (x,y) , i.e.,
if you are given a point (s,t) in isoparametric coordinates, then
you can compute the coordinates of the point in the (x,y)
coordinate system using the equations
x   N i ( s, t ) x i
i
y   N i ( s, t ) y i
i
2. The inverse map will never be explicitly computed.
8-noded Serendipity element
© 2002 Brooks/Cole Publishing / Thomson Learning™
t 7
4
1
1
1
6
8
1
1
3
5
2
s
8-noded Serendipity element: element shape functions in
isoparametric coordinates
1
1
N1 ( s, t )  (1  s )(1  t )(  s  t  1) N 5 ( s, t )  (1  t )(1  s)(1  s)
4
2
1
1
N 2 ( s, t )  (1  s )(1  t )( s  t  1) N 6 ( s, t )  (1  t )(1  t )(1  s)
4
2
1
1
N 3 ( s, t )  (1  s )(1  t )( s  t  1) N 7 ( s, t )  (1  t )(1  s)(1  s)
4
2
1
1
N 4 ( s, t )  (1  s )(1  t )(  s  t  1) N 8 ( s, t )  (1  s)(1  t )(1  t )
4
2
NOTES
1. Ni(s,t) is a simple polynomial in s and t. But Ni(x,y) is a complex
function of x and y.
2. The element edges can be curved in the mapped coordinates
3. A “midside” node in the parent element may not remain as a
midside node in the mapped element. An extreme example
t 5
2
1
1
1
8
6
1
3
7
y
1
4
s
2,6,3
1
5
8
7
x
4. Care must be taken to ensure UNIQUENESS of mapping
2
y
t
2
1
1
1
1
3
s
1
3
4
4
1
x
Isoparametric mapping in 2D: Triangular parent elements
t
2
1
s
3
Parent element: a right angled triangles
with arms of unit length
Key is to link the isoparametric
coordinates with the area coordinates
P(s,t)
t
1
1
1
2
t

2
s

2
1

(1  s  t )
2
s
 123 
 P 31
 P 23
 P12
L1 
 P 23
 s
 123
L2 
 P 31
t
 123
L3 
 P12
 1 s  t
 123
Now replace L1, L2, L3 in the formulas for the shape functions of
triangular elements to obtain the shape functions in terms of (s,t)
Example: 3-noded triangle
2 t
2 (x2,y2)
t
1
1
3
s
1
Parent shape functions
y
3 (x3,y3)
x
1 (x1,y1)
s
Isoparametric mapping
N1  s
x  N1 (s, t ) x1  N 2 (s, t ) x2  N 3 ( s, t ) x3
N2  t
y  N1 ( s, t ) y1  N 2 (s, t ) y 2  N 3 (s, t ) y 3
N3  1  s  t
Element matrices and vectors for a mapped 2D element
Recall: For each element
Displacement approximation
uNd
Strain approximation
ε Bd
Stress approximation
  DB d
k   e B D B dV
Element stiffness matrix
T
V
Element nodal load vector
f   e N X dV   e N T S dS
V
ST

 

 

T
f
STe
e
T
b
f
S
In isoparametric formulation
1. Shape functions first expressed in (s,t) coordinate system
i.e., Ni(s,t)
2. The isoparamtric mapping relates the (s,t) coordinates with the
global coordinates (x,y)
x   N i ( s, t ) x i
i
y   N i ( s, t ) y i
i
3. It is laborious to find the inverse map s(x,y) and t(x,y) and we do
not do that. Instead we compute the integrals in the domain of the
parent element.
NOTE
1. Ni(s,t) s are already available as simple polynomial functions
2. The first task is to find N i and N i
y
x
Use chain rule
N i ( x, y )
N i

s
x
N i ( x, y )
N i

t
x
N i
x

s
y
N i
x

t
y
y
s
y
t
In matrix form
 N i
 s
 N
i

 t

 x

 s
   x


 t

Can be
computed
 N i
 s
 N
i

 t
y   N i

s   x
y   N i

t  
 y
This is known as the
Jacobian matrix (J)
for the mapping
(s,t) → (x,y)
 N i 



 x 


J

 N 
i





 y 








We want to
compute these
for the B matrix
 N i

 x
 N
i


 y

 N i


1  s
  J  N
i



 t






How to compute the Jacobian matrix?
Start from
x  N ( s, t ) x

i
i
i
y   N i ( s, t ) y i
i
N ( s, t )
N ( s, t )
x
x
 i
xi ;
 i
xi
s
s
t
t
i
i
N i ( s, t )
N i ( s, t )
y
y

yi ;

yi
s
s
t
t
i
i
 x
 s
J 
x

 t
y 
s 
y 

t 
N i ( s, t )

xi

s
 i
 N i ( s, t ) x
i


t
 i
N i ( s, t )

yi 

s
i

N i ( s, t )
yi 


t
i

Need to ensure that det(J) > 0 for one-to-one mapping
3. Now we need to transform the integrals from (x,y) to (s,t)
Case 1. Volume integrals

1 1
Ve
f (s, t ) dV   e f (s, t ) h dA    f (s, t ) h det(J ) dsdt
A
1 1
h=thickness of element
This depends on the key result
dA  dxdy  det(J ) dsdt
Problem: Consider the following isoparamteric map
1
2
t
(6,6)
(3,6)
2
1
1
1
1
s
1
y
4
3
4
3
(3,1)
(6,1)
x
ISOPARAMETRIC COORDINATES
GLOBAL COORDINATES
Displacement interpolation
u  N1u1  N 2 u 2  N 3u3  N 4 u 4
v  N1v1  N 2 v2  N 3 v3  N 4 v4
Shape functions in isoparametric coord system
1
N 1 ( s, t )  (1  s )(1  t )
4
1
N 2 ( s, t )  (1  s )(1  t )
4
1
N 3 ( s, t )  (1  s )(1  t )
4
1
N 4 ( s, t )  (1  s )(1  t )
4
The isoparamtric map
x  N1 ( s, t ) x1  N 2 ( s, t ) x2  N 3 ( s, t ) x3  N 4 ( s, t ) x4
y  N1 ( s, t ) y1  N 2 ( s, t ) y2  N 3 ( s, t ) y3  N 4 ( s, t ) y4
(1  s )(1  t )
(1  s )(1  t )
(1  s )(1  t )
(1  s )(1  t )
6
3
3
6
4
4
4
4
3(1  s )

2
(1  s )(1  t )
(1  s )(1  t )
(1  s )(1  t )
(1  s )(1  t )
y
6
6
1
1
4
4
4
4
7  5t

2
x
In this case, we may compute the inverse map, but we will NOT do that!
The Jacobian matrix
y 
 x
 s
s 
J 


x

y



t 
 t

3

0
2



5
0


2


3(1  s)
2
7  5t
y
2
x
since
NOTE: The diagonal terms are due to stretching of the sides along the x-and ydirections. The off-diagonal terms are zero because the element does not shear.
2 / 3
J 1  
 0
0 
2 / 5

and
det( J ) 
15
4
Hence, if I were to compute the first column of the B matrix along the
positive x-direction
 N1

x


B1   0
 N
1

 y

I would use









 N1 
 N1 




2 / 3
 x 
1  s 

J

 N 


 0

N
1

1 




 t 


 y 

Hence
 N1   1  t 

 


x
6




 
B1   0    0 
 N  1  s 
1 



 y 
  10 

1 t  1 t 
 6 

 4 
 
0 




2 / 5  1  s  1  s 
 4 
 
 10 


The element stiffness matrix
k

V
B D B dV 
e
T
1

1
1 1
BT D B det( J )hdsdt