Transcript Finite Element Methods (FEM) Suzanne Vogel COMP 259 Spring, 2002 Definition of FEM The finite element method is the formulation of a global model to.

Finite Element Methods (FEM)
Suzanne Vogel
COMP 259
Spring, 2002
Definition of FEM
The finite element method is the formulation
of a global model to simulate static or dynamic
response to applied forces.
• Models: energy, force, volume,…
This differs from a mass spring system, which
is a local model.
Top-Down: Steps in FEM
1. Set up a global model in terms of the world
coordinates of mass points of the object.
These equations will be continuous.
2. Discretize the object into a nodal mesh.
3. Discretize the equations using finite
differences and summations (rather than
derivatives and integrals).
4. Use (2) and (3) to write the global equations
as a stiffness matrix times a vector of
(unknown) nodal values.
Top-Down: Steps in FEM
6. Solve for the nodal values.
•Static – nodal values at equilibrium
•Dynamic – nodal values at next time step
7. Interpolate values between nodal coordinates.
8
object
+
global model
discretize
interpolate
nodal mesh
+
local model
7
5
6
4 u 3
1
2
interpolate values between nodes
Bottom-Up: Steps in FEM
Nodes are point masses connected with springs.
A continuum equation is solved for the nodes,
and intermediate points are interpolated.
A collection of nodes forms an element.
A collection of elements forms the object.
8
5
7
6
4 u 3
1
2
Elements and Interpolations
Interpolating equations for an element are
determined by the number and distribution of
nodes within the element.
More nodes mean higher degree, for smoother
simulation.
Example: Hermite as 1D Cubic
Interpolation Equation
1. Assume
r (u)  au3  bu2  cu  d and
cubic equation
r (u)  N1 (u)  r1  N2 (u)  r2  N3 (u)  r3  N4 (u)  r4
equation using shape (blending) functions
r
u
Example: Hermite as 1D Cubic
Interpolation Equation
2. Normalize the element to [0,1] and rewrite
r (u)  au3  bu2  cu  d as a matrix equation
 r1  u
r  
 2   u
 r3  u
  
r4  u
3
1
3
2
3
3
3
4
2
1
2
2
2
3
2
4
u
u1
u
u2
u
u3
u
u4
or R0  U0  Q
0

1  a   1
   
1  b 

  27
1  c   8
    27
1 d 
 1
0
1
9
4
27
1
0
1
3
2
27
1
1
 a 
1  
b


1  c 
 d 
1
Example: Hermite as 1D Cubic
Interpolation Equation
3. Solve for the coefficients Q
1
R0  U0  Q  Q  U0  R0  M H R0
4. Plug the coefficients into the cubic equation
r (u)  au3  bu2  cu  d
5. Rewrite the cubic equation in the form
r (u)  N1 (u)  r1  N2 (u)  r2  N3 (u)  r3  N4 (u)  r4
Example: Hermite as 1D Cubic
Interpolation Equation
4 + 5. are equivalent to the steps
r (u )  U  Q  U  ( M H  R0 )  (U 0  M H )  R0
 r1 
r 
r (u )  N1 (u ) N 2 (u ) N 3 (u ) N 4 (u )  2 
 r3 
 
shape (blending) functions
r4 
values at the 4 nodes of the element
Example: Hermite as 1D Cubic
Interpolation Equation
r
shape (blending) functions within one element
Let Ni  Hi , u  t
u
U0
1
1D Elements
(x)
(x)
Example: bungee
(x)
2D Elements
(x,y)
(x,y)
Example: cloth
(x,y)
3D Elements
(x,y,z)
(x,y,z)
Example: skin
Static vs. Dynamic FEM
Static analysis is good for engineering, to find
just the end result.
Dynamic analysis is good for simulation, to
find all intermediate steps.
Types of Global Models[6]
Variational - Find the position function, w(t)
that minimizes the some variational integral.
This method is valid only if the position
computed satisfies the governing differential
equations.
Rayleigh-Ritz - Use the variational method
assuming some specific form of w(t) and
boundary conditions. Find the coefficients and
exponents of this assumed form of w(t).
Example of Variational Method[6]
Minimizing the variation w.r.t. w of the
variational function
 2


1
J ( w)     c2 w  c1 w w  2 fw d
2


under the conditions



c1  (   a)
 a w b w c
w
c2  (  b)
c3  (c  f )  0
satisfies the governing



 f (t )  0
equation, Lagrange’s  w  w
w
Equation
Types of Global Models[6]
Galerkin (weighted residual) - Minimize the
residual of the governing differential equation,
F(w,w’,w’’,…,t) = 0. The residual is the form of
F that results by plugging a specific form of the
position function w(t) into F. Find the
coefficients and exponents of this assumed form
of w(t).
Example of Galerkin Method[6]
We can approximate w(t) using Hooke’s Law
 1
1  L0

E  1
 L0
1
1
 
L0  f1 (t )  E   w1 (t ) 


1   f 2 (t )  E   w2 (t )
L0 
If we use that equation to compute the 1st and
2nd time derivatives of w, then we can compute
the residual as



 w  w 
 f (t )  0
w
Example of Static, Elastic FEM
Problem: If you apply the pressure shown, what
is the resulting change in length?
Object
Example of Static, Elastic FEM
First step. Set up a continuum model:
Infinitessimal length:
•F = force
P
 L 
•P = pressure
F   E

A
 L 
•A = area
•L = initial length Entire length:
•E = Young’s modulus
PL
L 
AE
PL
PL un 1
PL un 1
L  
du 
du 
du


u 0 AE
E u 0 wh(u )
wE u 0 h(u )
un
Example of Static, Elastic FEM
Since the shape is regular, we can integrate to
find the solution analytically. But suppose we
want to find the solution numerically.
Next step. Discretize the object.
Example of Static, Elastic FEM
n1
Discretization of object into
linear elements bounded by nodes
2
3
1
4
n2
n3
n4
n5
Example of Static, Elastic FEM
Next step. Set up a local model.
Stress-Strain Relationship (like Hooke’s Law)
Young’s modulus
i 
Ei , j
L0
distance between adjacent nodes
L  k i , j L  ki , j (ri  r j  L0 )  k i , j (ri  r j )  k i , j L0
stress (elastic force)
 j   i  ki, j (ri  rj )  ki, j L0
Example of Static, Elastic FEM
Next step. Set up a local (element) stiffness matrix.
 i  ki, j (ri  rj )  ki, j L0
 j  ki, j (ri  rj )  ki, j L0
nodal stresses
Rewrite the
above as a
matrix equation.
element stiffness matrix
nodal coordinates
 i   k i , j  k i , j   ri   k i , j L0 
   
    k


k
r
k
L
i, j   j 
 j   i, j
 i, j 0 
 i  k i , j L0   k i , j  k i , j   ri 


 


 j  k i , j L0   k i , j k i , j  r j 
Same for the
  j  k j , j 1 L0   k j , j 1



adjacent element.  j 1  k j , j 1 L0   k j , j 1
 k j , j 1   r j 
 

k j , j 1  r j 1 
Example of Static, Elastic FEM
Now, all of the element stiffness matrices are as
r is the x-coordinate of node u
follows.
i
i
 1  k1, 2 L0   k1, 2
  k L    k
 2 1, 2 0   1, 2
 k1, 2   r1   2  k 2,3 L0   k 2,3


k1, 2  r2   3  k 2,3 L0   k 2,3
 k 2,3  r2 
 

k 2,3   r3 
 3  k 3, 4 L0   k 3, 4
  k L    k
3, 4 0 
 4
 3, 4
 k 3, 4   r3   4  k 4,5 L0   k 4,5

  


k 3, 4  r4   5  k 4,5 L0   k 4,5
 k 4,5  r4 
 

k 4,5  r5 
n1
1
n2
2
n3
3
n4
4
n5
Example of Static, Elastic FEM
Next step. Set up a global stiffness matrix.
Pad the element stiffness matrices with zeros
and sum them up. Example:


 2
 3



0
 0
 k 2,3 L0  0 k 2,3
 k 2,3 L0   0  k 2,3
 
0
0
 0
 0
0
0
0
0
 k 2,3
k 2,3
0
0
0 0  r1 
0 0 r2 
0 0   r3 
  
0 0 r4 
0 0  r5 
Example of Static, Elastic FEM
Final step. Solve the matrix equation for the
Nodal coordinates.
nodal coordinates.
Solve for these!
Applied forces
 1  k1, 2 L0

  k1, 2
  (k  k ) L   k
1, 2
2,3
0
 2
 1, 2
 3  (k 2,3  k 3, 4 ) L0    0

 


(
k

k
)
L
3, 4
4,5
0
 4
 0

  0
 5  k 4,5 L0

 
 k1, 2
0
0
k1, 2  k 2,3
 k 2,3
 k 2,3
k 2 , 3  k 3, 4
0
 k 3, 4
0
 k 3, 4
k 4 , 5  k 3, 4
0
0
 k 4,5
Global stiffness matrix.
Captures material properties.
0   r1 
0  r2 
 
0    r3 
  
 k 4,5  r4 
k 4,5   r5 
Elastic FEM
A material is elastic if its behavior depends only
on its state during the previous time step.
•Think: Finite state machine
The conditions under which an “elastic”
material behaves elastically are:
•Force is small.
•Force is applied slowly and steadily.
Inelastic FEM
A material is inelastic if its behavior depends on
all of its previous states.
A material may behave inelastically if:
•Force is large - fracture, plasticity.
•Force is applied suddenly and released, i.e., is
transient - viscoelasticity.
Conditions for elastic vs. inelastic depend on
the material.
Examples of Elasticity
Elasticity
•Springs, rubber, elastic, with small, slowlyapplied forces
Examples of Inelasticity
Inelasticity
•Viscoelasticity
•Silly putty bounces under transient force (but
flows like fluid under steady force)
•Plasticity
•Taffy pulls apart much more easily under
more force (material prop.)
•Fracture
•Lever fractures under heavy load
Linear and Nonlinear FEM
Similarly to elasticity vs. inelasticity, there are
conditions for linear vs. nonlinear deformation.
Often these coincide, as in elastoplastic.

L
  Ee : e 
L0
=e
Elastic vs. Inelastic FEM
Elastic Deformation

e
loading
unloading
e
t

L
Hooke’s Law
or   Ee : e 
f  a0e
L
0
Young’s modulus
stress
strain
•Describes spring without damping
•Linear range of preceding stress vs. strain graph
Elastic vs. Inelastic FEM
Damped Elastic Deformation

e
loading
unloading
.
a1e
.
a1e
e
t
Rate of deformation is constant.

f  a1 e a0e
viscous
linear stress
Elastic vs. Inelastic FEM
Viscoelastic Deformation

depends on time t
loading
e
unloading
.
e
This graph is actually viscous,
but viscoelastic is probably similar


t
Rate of deformation is greatest
immediately after starting
loading or unloading.


b2f  b1 f  b0 
f  a2 e  a1 e a0e
new term! viscous linear stress
Elastic vs. Inelastic FEM
Elastoplastic Deformation
depends on force f

e
This graph is actually plastic,
but viscoelastic is probably similar

f
x
loading
x
loading
compare
f  a0e
e
unloading
e
x
Elastic vs. Inelastic FEM
Fracture
depends on force f
•Force response is locally discontinuous
•Fracture will propogate if energy release rate
is greater than a threshold
x

loading
e
x
unloading
Elastic vs. Inelastic
FEM4,5
1. World
coordinates w
in inertial frame

(a frame with
constant velocity)
2. Object
(material)
coordinates r
in non-inertial 
frame
r(w,t) = rref(w,t) + e(w,t)
ref
r
object, or
non-inertial
frame


world, or
inertial frame

origin of
= center of mass in

Elastic vs. Inelastic
FEM4,5
Transform
•reference
component rref
•elastic component e
•object frame 
w.r.t. world frame 
r(w,t) = rref(w,t) + e(w,t)
ref
r


Elastic vs. Inelastic FEM
All these equations are specific for:
•Elasticity
•Viscosity
•Viscoelasticity
•Plasticity
•Elastoplasticity
•Fracture
•(not mentioned) “Elastoviscoplasticity”
Ideally: We want a general equation that will
fit all these cases.
Elastic vs. Inelastic
FEM4,5
A More General Approach
To simulate dynamics we can use Lagrange’s
equation of strain force. At each timestep, the
force is calculated and used to update the
object’s state (including deformation).
Lagrange’s Equation
elastic potential energy

f ( w, t )   w  w
w


mass density damping density
stress component
of force
 ( / L0 ) E  e 

  
w
w
L0  w 
matrices
Elastic vs. Inelastic
Given:
Mass density and damping density FEM
are known.
4,5
Elastic potential energy derivative w.r.t. r can be
approximated using one of various equations.
vector
vector
next slide
The current position wt of all nodes of the object
are known.
Unknown:
The new position wt+dt of nodes is solved for at
each timestep.


Lagrange’s

f ( w, t )   w  w
Equation
w
Elastic vs. Inelastic
FEM
4,5



elastic potential energy
f ( w, t )   w  w
w
For both elastic and inelastic deformation,
express elastic potential energy as an integral
in terms of elastic potential energy density.
elastic potential energy density
 e  2e  3e 
     r , , 2 , 3 ,...dr

 r r r

Elastic vs. Inelastic
FEM
Elastic potential energy density  can be 4,5
approximated using one of various equations
which incorporate material properties.
• Elastic deformation: Use tensors called metric
(1D, 2D, 3D stretch), curvature (1D, 2D
bend), and “twist” (1D twist).
• Inelastic deformation: Use controlledcontinuity splines.
Elastic FEM4
For elastic potential energy density in 2D, use
• metric tensors G (for stretch)
• curvature tensors B (for bend)
 (r) || G  G ||  || B  B || 
0
2
0
2
|| M || = weighted norm of matrix M
Elastic FEM4
Overview of derivation of metric tensor
Since the metric tensor G represents stretch, it
incorporates distances between adjacent points.
 w   w 
dri drj
  
dL  dw dw   


i , j 1, 2  ri   r j 
2
object coordinates

G
i , j 1, 2
i, j
dri drj  G1,1
 G1,1
 1 1  
G2,1
 dr1dr1 
 dr dr 
G2, 2   1 2 
 dr2 dr1 


dr
dr
 2 2
world coordinates
G1, 2
G2,1
G1, 2   dr1   dr1 




G2, 2  dr2  dr2 
T
Elastic FEM4
Overview of metric and curvature tensors.
From the previous slides, we found:
 w   w 
 represents stretch
Gi , j ( w(r ))     
 r 

r
 i  j
Similarly:
 2w
Bi , j ( w(r ))  
 r r
 i j




represents bend
Theorem. G and B together determine shape.
Elastic FEM4
For elastic FEM, elastic potential energy
density in 2D incorporates changes in the
metric tensor G and the curvature tensor B.
 (r) || G  G ||  || B  B || 
0
2
0
2
weights = material properties
|| M || = weighted norm of matrix M
Inelastic FEM5
For inelastic FEM, elastic potential energy
density is represented as a controlledcontinuity spline. weighting function = material property


1
m!
  
w j  j1 j 2
2 m0 | j|m j1! j 2 !...j d !  r r ...r jd
p
m
j

e


2
For some degree p, dimensionality d, compute
the sum of sums of all combinations of
weighted 1st, 2nd,…, mth derivatives of strain e
w.r.t. node location r, where m <= p.
Inelastic FEM5
Then the elastic potential energy density
derivative w.r.t. strain e is:
weighting function = material property
 mj

   1m  j1 j 2
jd
e m0
| j|  m r r ...r
p
m



m!
j


w
 j1! j 2 !...j d ! j  r j1r j 2 ...r jd



e


Example: p = 2, d = 3

0!


w00 e 
e 0!0!
r1
2

r1 r2
 1!
   



 1!0! w10  r e    r
2
 1 

 2!
 2
  2
 w11 

 1!1!  r r e    r 2
1
 1 2 

 1!
  



 0!1! w01  r e  
 2 

 2!
 2


 2!0! w20  r 2
 1

  2
e    2
  r2
 2!
 2


 0!2! w02  r 2
 2


e  

Recap
Elastic vs. Inelastic
 Lagrange’s Eq’n
FEM
f ( w, t )   w  w
4,5
w

total force
(includes stress)

elastic
potential energy
Elastic

How it has been
expanded and is continuing
to be expanded...
 e  2 e  2 e 
     r , , 2 , 2 ,...dr
u
 r r r

4
 (r) || G  G 0 ||2  || B  B 0 || 2
5
material properties
elastic potential
energy density
Inelastic


1
m!

  
wj
2 m0 | j|m j1! j 2 !...j d !  r j1r j 2 ...r jd
p
m
j
5
p
j

m
   1  j1 j 2
jd
e m0

r

r
...

r
| j|  m
m
m



m
!
j

w j  j1 j 2
 j1! j 2 !...j d !  r r ...r jd



e


2

e


Elastic FEM4
Continuing

f ( w, t )   w  w
w

 e  2 e  2 e 
     r , , 2 , 2 ,...dr
u
 r r r

elastic
potential
energy
 (r )  


 (r) || G  G 0 ||2  || B  B 0 || 2
  G  G 
0 2
i , j 1, 2
i, j
  i , j B  B
 dr
0 2
i, j (r, w)  i, j (w)Gi, j  Gi0, j 
i, j (r, w)   i, j (w)Gi, j  Gi0, j 
>0: surface wants to shrink
<0: surface wants to expand
>0: surface wants to flatten
<0: surface wants to bend
Inelastic FEM5
Continuing
elastic potential
energy density

f ( w, t )   w  w
w

elastic potential energy
 e  2 e  2 e 
     r , , 2 , 2 ,...dr
u
 r r r

 mj

   1m  j1 j 2
jd
e m0
| j|  m r r ...r
p

m



m!
j


w
 j1! j 2 !...j d ! j  r j1r j 2 ...r jd


strain
Deformation has been modeled by
approximating elastic potential energy.

e


Inelastic FEM5
Continuing
Now rigid-body motion and other aspects of
deformation must be computed using physics
equations of motion.
In this way, both (in)elastic deformation and
rigid-body motion can be modeled, providing a
very general framework.
r(w,t) = rref(w,t) + e(w,t)
Inelastic FEM5
Motion of object (non-inertial) frame w.r.t.
world (inertial) frame

c(t )    (r ) w(r , t ).dr


f ( w, t )   w  w
w


d
d
f  (m c) 
dt
trans dt

d
d
f  (I  ) 
dt
rot dt






w(r , t )  c(t )   (t )  c(t )  e(r , t )

v




 e(r , t ).dr    w(r , t ).dr



r  e .dr   r  w .dr

Combines
dynamics of
deformable
and rigid
bodies








d

f   (t )  (  e)   c     (  r )  2   e    r   w
dt
e
elastic
e

f ( w, t )   w  w
w


Inelastic FEM5
Velocity of node of object (non-inertial) frame
w.r.t. world (inertial) frame (radians / sec) x (radius)




w(r , t )  c(t )   (t )  c(t )  e(r , t )
w.r.t. world
velocity of reference velocity of elastic
component
component
c(t )
e(t )

w(r , t )
 (t )

Identically, in another
coordinate system,
r(w,t) = rref(w,t) + e(w,t)
w.r.t. object

f ( w, t )   w  w
w

f
rot


d
d
 (I  ) 
dt
dt
Inelastic FEM5




r  e .dr   r  w .dr

angular momentum
  ( w22  w32 )
w1 w2 )
 w1 w3 


2
2
I (t )    w2 w1
 ( w1  w3 )  w2 w3 dw(t )

  w w
2
2 
 w3 w2
 ( w2  w3 ) 
3 1
inertia tensor
Angular momentum is conserved in the absense
of force. So a time-varying angular momentum
indicates the presence of foce.

f ( w, t )   w  w
w

f
rot


d
d
 (I  ) 
dt
dt
Inelastic FEM5




r  e .dr   r  w .dr

indicates changing angle between position and direction of stretch
r (t )

e(t )

f ( w, t )   w  w
w


Inelastic FEM5



 

 
d 
f   (t )  (  e)   c     (  r )  2   e    r   w
dt
e
elastic
e
inertial
centripetal
Coriolis
transverse
damping
restoring
elastic potential energy
strain
If the reference component has no translation or

rotation, then e
d

f   (t ) 
dt
(  e) 
e
Furthermore, if the elastic component has no
acceleration, then f e   (t )  
e

f ( w, t )   w  w
w


Inelastic FEM5
Recall that non-elastic behavior is characterized
by acceleration of the elastic component
(strain)...

d

e
f   (t )  (  e) 
dt
e
And elastic behavior is characterized by
constant velocity of strain.

e
f   (t ) 
e

x
loading
f  a0e
e
Elastic vs. Inelastic
4,5
Now Lagrange’s equation has beenFEM
expanded.

f ( w, t )   w  w
w


Final Steps
•Discretize using finite differences (rather than
derivatives).
•Write as a matrix times a vector of nodal
coordinates (rather than a single mass point).
•Solve for the object’s new set of positions of
all nodes.
Discretization of FEM4,5
Discretize Lagrange’s equation over all nodes

f ( w, t )   w  w
w



f (t )  M  w C  w
w


Procedure described in [4] but not [5]
Discretization of Elastic FEM4
2w
w
M 2 C
 K ( wt  t ) wt  t  f t  t
t
t
 MDt ( Dt ( wt ))  CDt ( wt )  K ( wt ) wt  t
wt  t  2wt  wt  t
wt  t  wt  t
M
C
 K ( wt ) wt  t
2
2t
t
1
1 
 1

 2

 1
 2 M 
 K ( wt )  wt  t  f t   2 M  wt   2 M 
C  wt  t
2t
2t 
 t

 t

 t
1 
1  wt  wt  t
 3
1
 ft   2 M 
C  wt   M  C 
2t 
2 
t
 t
 t
1 
1 
 3
1
 ft   2 M 
C  wt   M  C vt
2t 
2 
 t
 t
 At wt  t  g t ( wt , wt  t ) _ where
1
1 
1 
 1

 3
1
At   2 M 
 K ( wt ) , g t ( wt , wt  t )   2 M 
C  wt   M  C vt
2t
2t 
2 
 t

 t
 t
Results of Elastic FEM4
Results of Elastic FEM4
Results of Elastic FEM4
Results of Inelastic FEM5
3D plasticine bust of Victor Hugo.
180 x 127 mesh; 68,580 equations.
Results of Inelastic FEM5
Sphere pushing through 2D mesh.
23 x 23 mesh; 1,587 equations.
Yield limit is uniform, causing linear tears.
Results of Inelastic FEM5
2D paper tearing by opposing forces.
30 x 30 mesh; 2,700 equations.
Yield limit is perturbed stochastically,
causing randomly-propogating tears.
References
0. David Baraff. Rigid Body Simulation.
Physically Based Modeling, SIGGRAPH
Course Notes, August 2001.
1. George Buchanan. Schaum’s Outlines:
Finite Element Analysis. McGraw-Hill, 1995.
2. Peter Hunter and Andrew Pullan. FEM/BEM
Notes. The University of Auckland, New
Zealand, February 21 2001.
References
3. Tom Lassanske. [Slides from class lecture]
4. Demetri Terzopoulost, John Platt, Alan Barr,
and Kurt Fleischert. Elastically Deformable
Models. Computer Graphics, Volume 21,
Number 4, July 1987.
5. Demetri Terzopoulos and Kurrt Fleiseher.
Modeling Inelastic Deformation:
Viscoelasticity, Plasticity, Fracture. Computer
Graphics, Volume 22, Number 4, August 1988
Notation
r  object_ coordinates
w  world _ coordinates
e  strain _(stretch)
  stress _( force)
E  Young' s _ modulus
  elastic _ potential_ energy
  elastic _ potential_ energy_ density