Transcript MANE 4240 & CIVL 4240 Introduction to Finite Elements Prof. Suvranu De FEM Discretization of 2D Elasticity.

MANE 4240 & CIVL 4240
Introduction to Finite Elements
Prof. Suvranu De
FEM Discretization of
2D Elasticity
Reading assignment:
Lecture notes
Summary:
• FEM Formulation of 2D elasticity (plane stress/strain)
•Displacement approximation
•Strain and stress approximation
•Derivation of element stiffness matrix and nodal load vector
•Assembling the global stiffness matrix
• Application of boundary conditions
• Physical interpretation of the stiffness matrix
Recap: 2D Elasticity
Volume
element dV
Xb dV
Xa dV
v
y
py
px
Volume (V)
u
ST
x
Su
x
Examples: concept of displacement field
Su: Portion of the
boundary on which
displacements are
prescribed (zero or
nonzero)
ST: Portion of the
boundary on which
tractions are prescribed
(zero or nonzero)
Example
y
2
1
2
2
4
3
x
For the square block shown above, determine u and v for the
following displacements
Case 2: Pure shear
y
y
Case 1: Stretch
1/2
2
4
2
1
x
Solution
Case 1: Stretch
ux
y
v
2
Check that the new coordinates (in the deformed configuration)
x'  x  u  2x
y
'
y  yv 
2
Case 2: Pure shear
u  y/4
v0
Check that the new coordinates (in the deformed configuration)
x'  x  u  x  y / 4
y'  y  v  y
Recap: 2D Elasticity
Displacement field u  u ( x, y)
Strain - Displacement Relation  u
Stress - Strain Law   D  Du
u (x, y)
u

v (x, y)


0

x


x 
 x 




0
 
 

   y 
   y 
y 

 
 



 xy 
 xy 
 y x 
For plane stress
For plane strain
(3 nonzero stress components) (3 nonzero strain components)




1



0
1 
0 



E
E
 
D
1 
0 


D

1
0
1   1  2  
1  2 
1  2 
1  
0
0


0
0


2


2 

Strong formulation
Equilibrium equations
   X  0 in V
T
Boundary conditions
1. Displacement boundary conditions: Displacements are specified on
portion Su of the boundary
u u
specified
on Su
2. Traction (force) boundary conditions: Tractions are specified on
portion ST of the boundary
Now, how do I express this mathematically?
But in finite element analysis we DO NOT work with the strong
formulation (why?), instead we use an equivalent Principle of
Minimum Potential Energy
Principle of Minimum Potential Energy (2D)
Definition: For a linear elastic body subjected to body forces
X=[Xa,Xb]T and surface tractions TS=[px,py]T, causing
displacements u=[u,v]T and strains  and stresses , the potential
energy P is defined as the strain energy minus the potential energy
of the loads (X and TS)
PU-W
Volume
element dV
Xb dV
Xa dV
v
y
Volume (V)
u
ST
x
Su
x
U
1
T

 dV

2 V
W   u X dV   u T S dS
T
V
T
ST
py
px
Strain energy of the elastic body
Using the stress-strain law   D 
1
1
T
T
U     dV    D  dV
2 V
2 V
In 2D plane stress/plane strain
1
T
U     dV
2 V
T
 x    x 
1    
   y    y  dV
2 V   
 xy   xy 
1
   x x   y y   xy xy  dV
2 V
Principle of minimum potential energy: Among all admissible
displacement fields the one that satisfies the equilibrium equations
also render the potential energy P a minimum.
“admissible displacement field”:
1. first derivative of the displacement components exist
2. satisfies the boundary conditions on Su
Finite element formulation for 2D:
Step 1: Divide the body into finite elements connected to each
other through special points (“nodes”)
py
v3
3
px
4
3
u3
u 1 
v4
2
v 
v2
Element ‘e’
v
1 

1
4
u 2 
u
u4
ST
 
v1
2 u2
v 2 
y
d 
x
y
u3 

Su
u1
v 3 
1
 
v
x
x
u 4 
u
v 
 4
Total potential energy
1
T
T
T
P     dV   u X dV   u T S dS
V
ST
2 V
Potential energy of element ‘e’:
1
T
T
T
P e     dV   u X dV   u T S dS
V
S
2 V
e
e
e
T
This term may or may not be present
depending on whether the element is
actually on ST
Total potential energy = sum of potential energies of the elements
P  Pe
e
Step 2: Describe the behavior of each element (i.e., derive the
stiffness matrix of each element and the nodal load vector).
Inside the element ‘e’
v3
(x3,y3) 3
v4
v2
Displacement at any point x=(x,y)
u (x, y)
u

v (x, y)
(x4,y4)
Nodal displacement vector
4
u4
u 1 
u2
v
v1
2
v  where
1 
u
(x
,y
)

u1=u(x1,y1)
2
2
y
u 2  v =v(x ,y )
1
1 1
u


1
1(x1,y1)
v 2  etc
d 
x
u 3 
v 3 
 
u 4 
v 
 4
u3
Recall
Strain - Displacement Relation  u
Stress - Strain Law   D  Du
 x 
 
   y 
 
 xy 
 x 
 
   y 
 
 xy 


 x
0



 y

0


y 


x 
If we knew u then we could compute the strains and stresses within the
element. But I DO NOT KNOW u!!
Hence we need to approximate u first (using shape functions) and
then obtain the approximations for  and  (recall the case of a 1D bar)
This is accomplished in the following 3 Tasks in the next slide
TASK 1: APPROXIMATE THE DISPLACEMENTS WITHIN
EACH ELEMENT
Displacement approximation in terms of shape functions
uNd
TASK 2: APPROXIMATE THE STRAIN and STRESS WITHIN
EACH ELEMENT
Strain approximation
ε Bd
Stress approximation
  DB d
TASK 3: DERIVE THE STIFFNESS MATRIX OF EACH
ELEMENT USING THE PRINCIPLE OF MIN. POT ENERGY
We’ll see these for a generic element in 2D today and then derive
expressions for specific finite elements in the next few classes
TASK 1: APPROXIMATE THE DISPLACEMENTS
WITHIN EACH ELEMENT
Displacement approximation in terms of shape functions
v3
3
u3
v4
v2 Displacement approximation within element ‘e’
4
u4
u2
v1
2
u (x,y)  N1(x,y) u1  N 2 (x,y) u 2  N 3(x,y) u 3  N 4 (x,y) u 4
y
u1
v (x,y)  N1(x,y) v1  N 2 (x,y) v2  N 3(x,y) v3  N 4 (x,y) v4
1
v
x
u
u (x,y)  N1(x,y) u1  N 2 (x,y) u 2  N 3(x,y) u 3  N 4 (x,y) u 4
v (x,y)  N1(x,y) v1  N 2 (x,y) v2  N 3(x,y) v3  N 4 (x,y) v4
u (x, y)  N 1
u

v (x, y)  0
uNd
0
N2
0
N3
0
N4
N1
0
N2
0
N3
0
u 1 
v 
 1
u 2 
 
0  v 2 
 

N 4  u 3 
v 3 
 
u 4 
v 
 4
We’ll derive specific expressions of the shape functions for
different finite elements later
TASK 2: APPROXIMATE THE STRAIN and STRESS WITHIN
EACH ELEMENT
Approximation of the strain in element ‘e’
N 3 (x, y)
N 2 (x, y)
N 4 (x, y)
u (x, y) N1(x, y)
x 

u1 
u2 
u3 
u4
x
x
x
x
x
N 3 (x, y)
N 2 (x, y)
N 4 (x, y)
v (x, y) N1(x, y)
y 

v1 
v2 
v3 
v4
y
y
y
y
y
N1(x, y)
u (x, y) v (x, y) N1(x, y)
 xy 


u1 
v1  ......
y
x
y
x
x 
 
   y 
 xy 
 
u 1 
v 
 N 1(x, y)
 1 
N 3 (x, y)
N 2 (x, y)
N 4 (x, y)
0
0
0
0

 u 2 

x

x

x

x

 

N
(x,
y)

N
(x,
y)

N
(x,
y)

N
(x,
y)
3
1
2
4
 v 2 

0
0
0
0
 

y
y
y
y  u 3 
 N (x, y) N (x, y) N (x, y) N (x, y) N (x, y) N (x, y) N (x, y) N (x, y)
3
3
1
2
2
4
4
 1
 v 3 
x
y
x
y
x
y
x  u 
 y
  4 
B
v 
 4
ε Bd
Compact approach to derive the B matrix:
Displacement field u  N d
Strain - Displacement Relation  u   N d  Bd
B N
Stress approximation within the element ‘e’
Stress - Strain Law   D
  D B 
TASK 3: DERIVE THE STIFFNESS MATRIX OF EACH
ELEMENT USING THE PRINCIPLE OF MININUM
POTENTIAL ENERGY
Potential energy of element ‘e’:
1
T
T
T
P e   e   dV   e u X dV   e u T S dS
V
ST
2 V
Lets plug in the approximations
uNd
ε Bd
  DB d
1
T
T
T
P e (d)   e DB d  B d  dV   e N d  X dV   e N d  T S dS
V
ST
2 V
Rearranging
1 T
T
T
T
T
T
P e (d)  d   e B D B dV  d  d   e N X dV   d   e N T S dS 

V

 ST

2 V
1 T
T
T
T
T
 d   e B D B dV  d  d   e N X dV   e N T S dS 
ST
V
2 V

k
f
1 T
T
d k dd f
2
From the Principle of Minimum Potential Energy
 P e (d ) 
P e (d)
k d f 0
d
Discrete equilibrium equation for element ‘e’
kd f
Element stiffness matrix for element ‘e’
k   e B D B dV
T
V
For a 2D element, the size of the k matrix is
2 x number of nodes of the element
Question: If there are ‘n’ nodes per element, then what is the size of
the stiffness matrix of that element?
Element nodal load vector
f   e N X dV   e N T S dS
V
ST

 

 

T
f
T
b
Due to body force
f
S
Due to surface traction
STe
e
If the element is of thickness ‘t’
k   e t B D B dA
T
A
For a 2D element, the size of the k matrix is
2 x number of nodes of the element
dA
dV=tdA
Element nodal load vector
t
f   e t N X dA   e t N T S dl
A
lT



 



T
T
f
f
b
Due to body force
S
Due to surface traction
The properties of the element stiffness matrix
k   e B D B dV
T
V
1. The element stiffness matrix is singular and is therefore noninvertible
2. The stiffness matrix is symmetric
3. Sum of any row (or column) of the stiffness matrix is zero!
(why?)
Computation of the terms in the stiffness matrix of 2D elements
v4
v3
4
The B-matrix (strain-displacement) corresponding to this element is
3
u4
u3
u1
v
y
v1 (x,y)
v2
u
2
1 u1
u2
 N1 (x,y)

x


0


 N1 (x,y)

y

v1
u2
v2
0
N 2 (x,y)
x
0
N1 (x,y)
y
N1 (x,y)
x
0
N 2 (x,y)
y
N 2 (x,y)
y
N 2 (x,y)
x
u3
N 3 (x,y)
x
0
N 3 (x,y)
y
x
We will denote the columns of the B-matrix as


 N1 (x,y) 


0




x


 N1 (x,y)  ; and so on...
B u1  
0
 ; B v1  

y
 N (x,y) 


1



N
(x,y)


1
y




x
v3
0
N 3 (x,y)
y
N 3 (x,y)
x
u4
v4
N 4 (x,y)
x



N 4 (x,y) 

y

N 4 (x,y) 

x

0
N 4 (x,y)
y
0
The stiffness matrix corresponding to this element is
k   e B D B dV
T
which has the following form
V
u1
v1
 k11
k
 21
 k31

k
k   41
 k51

 k61
k
 71

 k81
k12
k 22
k32
k 42
k52
k62
k72
k82
u2
k13
k 23
k33
k 43
k53
k63
k73
k83
v2
u3
k14
k 24
k34
k 44
k54
k64
k74
k84
k15
k 25
k35
k 45
k55
k65
k75
k85
v3
k16
k 26
k36
k 46
k56
k66
k76
k86
u4
v4
k17
k 27
k 37
k 47
k57
k 67
k 77
k87
k18 
k 28 

k 38 

k 48 
k58 

k 68 
k 78 

k88 

u1
v1
u2
v2
u3
v3
u4
v4
The individual entries of the stiffness matrix may be computed as follows
k11   e Bu1 D Bu1 dV; k12   e Bu1 D Bv1 dV; k13   e Bu1 D Bu2 dV,...
T
T
V
V
V
k21   e Bv1 D Bu1 dV; k21   e Bv1 D Bv1 dV;.....
T
V
T
V
T
Step 3: Assemble the element stiffness matrices into the global
stiffness matrix of the entire structure
For this create a node-element connectivity chart exactly as in 1D
v3
Element #1
3
v1
v4
1
v2
u1
y
u3 ELEMENT Node 1 Node 2 Node 3
4
u4
Element #2
u2
2
v
u
x
1
1
2
3
2
2
3
4
Stiffness matrix of element 1
u1 v1 u2 v2 u3 v3
k
(1)


















Stiffness matrix of element 2
u2 v2 u3 v3 u4 v4


u1




v1


( 2)
u2 k  



v2


u3




v3
There are 6 degrees of freedom (dof) per element (2 per node)
u2
v2
u3
v3
u4
v4
k
(1)
Global stiffness matrix
u1 v1 u2 v2 u3 v3 u4 v4




K  




 u1
 v1
 u
 2
v2

 u3
 v3
 u4
 8v48
How do you incorporate boundary conditions?
Exactly as in 1D
k
( 2)
Finally, solve the system equations taking care of the
displacement boundary conditions.
Physical interpretation of the stiffness matrix
Consider a single triangular element. The six corresponding
equilibrium equations ( 2 equilibrium equations in the x- and ydirections at each node times the number of nodes) can be written
symbolically as
kd f
v1
1
u1
v2
v3
3
y
u2
2
x
u3
 k11
k
 21
k 31

k 41
k 51

k 61
k12
k 22
k 32
k 42
k 52
k 62
k13
k 23
k 33
k 43
k 53
k 63
k14
k 24
k 34
k 44
k 54
k 64
k15
k 25
k 35
k 45
k 55
k 65
k16   u 1  f1x 
 
k 26   v1  f1y 
k 36  u 2  f 2x 
    
k 46  v 2  f 2y 
k 56  u 3  f 3x 
   
k 66   v 3  f 3y 
Choose u1 = 1 and rest of the nodal displacements = 0
 k11  f1x 
k   f 
 21   1y 
k 31  f 2x 
  
k 41  f 2y 
k 51  f 3x 
   
k 61  f 3y 
u1=1
1
3
y
2
x
Hence, the first column of the stiffness matrix represents the
nodal loads when u1=1 and all other dofs are fixed. This is the
physical interpretation of the first column of the stiffness matrix.
Similar interpretations exist for the other columns
“Force” at d.o.f ‘i’ due to unit displacement at d.o.f ‘j’
k ij = keeping
all the other d.o.fs fixed
Now consider the ith row of the matrix equation k d  f
ki1
ki 2 ki3 ki 4 ki5 ki 6  d  f ix
This is the equation of equilibrium at the ith dof
Consistent and Lumped nodal loads
Recall that the nodal loads due to body forces and surface tractions
f b   e N X dV;
T
V
fS 
ST e
T
N T S dS
These are known as “consistent nodal loads”
1. They are derived in a consistent manner using the Principle of
Minimum Potential Energy
2. The same shape functions used in the computation of the stiffness
matrix are employed to compute these vectors
p per unit area
Example
y
1
b
2
b
3
x
Traction distribution on the 12-3 edge
px = p
py= 0
We’ll see later that
y(b  y)
b2  y 2
y(b  y)
N1 
; N2 
; N3  
2
2
2b
b
2b 2
N1
N2
N3
The consistent nodal loads are
y (b  y )
pb
F1x   p N1 dy   p
dy 
2
b
b
3
2b
b
b
b2  y2
4 pb
F2 x   p N 2 dy   p
dy 
2
b
b
3
b
b
b
y (b  y )
pb
F3 x   p N 3 dy    p
dy 
2
b
b
3
2b
b
b
y
b
1
pb/3
2
4pb/3
b
3
pb/3
x
The lumped nodal loads are
pb
F1x 
2
F2 x  pb
pb
F3 x 
2
y
b
1
pb/2
2
pb
x
b
3
pb/2
Lumping produces poor results and will not be pursued further
Summary: For each element
Displacement approximation in terms of shape functions
uNd
Strain approximation in terms of strain-displacement matrix
ε Bd
Stress approximation
  DB d
Element stiffness matrix
k   e B D B dV
T
V
Element nodal load vector
f   e N X dV   e N T S dS
V
ST

 

 

T
f
T
b
f
S