Transcript Finite Element Method

The Finite Element Method
A Practical Course
CHAPTER1:
FOUNDMENTALS FOR FINITE ELEMENT
METHOD
CONTENTS



STRONG AND WEAK FORMS OF GOVERNING EQUATIONS
HAMILTON’S PRINCIPLE
FEM PROCEDURE
–
–
–
–
–
–
–




Domain discretization
Displacement interpolation
Formation of FE equation in local coordinate system
Coordinate transformation
Assembly of FE equations
Imposition of displacement constraints
Solving the FE equations
STATIC ANALYSIS
EIGENVALUE ANALYSIS
TRANSIENT ANALYSIS (reading materials)
REMARKS
STRONG AND WEAK FORMS OF
GOVERNING EQUATIONS

System equations: strong form (PDE), difficult to
solve.
 Weak (integral) form: requires weaker continuity
on the dependent variables (e.g., u, v, w).
 Weak form is often preferred for obtaining an
approximated solution.
 Formulation based on a weak form leads to a set
of algebraic system equations – FEM.
HAMILTON’S PRINCIPLE

“Of all the admissible time histories of
displacement the most accurate solution makes the
Lagrangian functional a minimum.”

An admissible displacement must satisfy:
– The compatibility conditions
– The essential or the kinematic boundary conditions
– The conditions at initial (t1) and final time (t2)
HAMILTON’S PRINCIPLE

Lagrangian functional
Mathematically
  t Ldt  0 where L=T-P+Wf
t2
1
1
T   U TUdV
2V
1
1
T
Π   ε σ dV   ε T cε dV
2V
2V
W f   U T f b dV 
V

Sf
U T f s dS f
(Kinetic energy)
(Potential energy)
(Work done by
external forces)
FEM PROCEDURE







Step 1: Domain discretization
Step 2: Displacement interpolation
Step 3: Formation of FE equation in local coordinates
Step 4: Coordinate transformation
Step 5: Assembly of FE equations
Step 6: Imposition of displacement constraints
Step 7: Solving the FE equations
Nodes
Triangular elements
Step 1: Domain discretization





The solid body is divided into Ne elements with proper
connectivity – compatibility.
All the elements form the entire domain of the problem
without any gap or overlapping – compatibility.
There can be different types of element with different
number of nodes.
The density of the mesh depends upon the accuracy
requirement of the analysis.
The mesh is usually not uniform, and a finer mesh is often
used in the area where the displacement gradient is larger.
Step 2: Displacement interpolation
y, v

Bases on local coordinate system,
the displacement within element is
interpolated using nodal displacements.
nd
3 (x3, y3)
(u3, v3)
fsy
fsx
A
2 (x2, y2)
(u2, v2)
1 (x1, y1)
(u1, v1)
x, u
U( x, y, z )   N i ( x, y, z ) di  N( x, y, z )d e
i 1
 d1
d
 2
di  

d n f

  displacement compenent 1

  displacement compenent 2


  displacement compenent n f

nf: Degree of freedoms at a node
 d1  displacements at node 1
d 
 2   displacements at node 2
de   
 
d nd   displacements at node nd
 
nd: number of nodes in an element
Step 2: Displacement interpolation

N is a matrix of shape functions
N( x, y, z )   N1 ( x, y, z ) N 2 ( x, y, z )
N nd ( x, y, z ) 

for node 1

for node nd
 N i1
 0
where N i  
 0

 0

for node 2
0
0
0
Ni2 0
0
0 
0
0
0 N in f






Shape function
for each
displacement
component at a
node
Displacement interpolation

Constructing shape functions
– Consider constructing shape function for
a single displacement component
– Approximate in the form
Basis function
u h ( x) 
nd

i 1
pi (x)i  pT (x)α
αT ={1 , 2 , 3 , ......, nd }
pT(x)={1, x, x2, x3, x4,..., xp}
(1D)
Pascal triangle of monomials: 2D
1 Constant terms: 1
3 terms
y
x
xy
x2
x
x
5
x 3y
4
x4y
x2y2
x3y2
y2
xy2
x 2y
x3
Linear terms: 2
6 terms
10 terms
Quadratic terms: 3
y3 Cubic terms: 4
xy3
x2y3
y4
xy4
Quartic terms: 5
y5
Quintic terms: 6
pT (x)  pT ( x, y )  1, x, y, xy, x 2 , y 2 ,..., x p , y p 
15 terms
21 terms
Pascal pyramid of monomials : 3D
1
Constant term: 1
4 terms
y
x
10 terms
Linear terms: 3
z
20 terms
y2
xy
x2
z2
2
xy2
xy
x3
xyz
xz2
x3y
x
x2yz
3
xz
2 2
xz
xz
3
y3
Cubic terms: 10
zy2
x2z
4
Quadratic terms: 6
yz
xz
35 terms
yz2
3
z
x2y2
xyz
xy3
2
xy2z
y4
zy3
Quartic terms: 15
2 2
zy
3
z4 z y
pT (x)  pT ( x, y, z )  1, x, y, z, xy, yz, zx, x 2 , y 2 , z 2 ,..., x p , y p , z p 
Displacement interpolation
– Enforce approximation to be equal to the nodal
displacements at the nodes
di = pT(xi)
i = 1, 2, 3, …,nd
or
 d1 
 p T (x1 ) 
de=P 
d 
 T

p
(
x
)

2 
2 
where
de =   ,
P 
 
 d nd 
 


 T

p (x nd ) 
Moment matrix
Displacement interpolation
– The coefficients in  can be found by
-1
α  P de
– Therefore, uh(x) = N( x) de

N(x)  pT (x)P -1  pT (x)P1-1 pT (x)P2-1

N2 ( x )
 N1 ( x )
  N1 (x) N 2 (x)
N n ( x) 

pT (x)Pn-1 

Nn ( x )

Displacement interpolation

Sufficient requirements for FEM shape
functions
1
1. Ni  x j    ij  
0
n
2.
 N ( x)  1
i 1
i
nd
3.
 N ( x) x
i 1
i
i
x
i  j,
j  1, 2,
i  j, i, j  1, 2,
, nd
, nd
(Delta function
property)
(Partition of unity property –
rigid body movement)
(Linear field reproduction property)
Step 3: Formation of FE equations in local
coordinates
Since U= Nde
Therefore, e = LU
Π
Strain matrix
 e = L N de= B de
1
1
1 T
T
T T
T


ε
c
ε
d
V
d
B
c
Bd
d
V
d
(
B
cBdV )d e

 e
e
e 
2 Ve
2 Ve
2
Ve
1 T
or P  d e k ed e where k   BT cBdV
e
2
Ve
(Stiffness matrix)
Step 3: Formation of FE equations in local
coordinates
Since U= Nde

U  Nd e
T
1
1
1 T
T
T T
T

U
U
d
V


d
N
Nd
d
V

d
(

N
NdV )de
e
e
e 


2 Ve
2 Ve
2
Ve
or
1 T 
T  d e m ed e
2
where
m e   NT NdV
Ve
(Mass matrix)
Step 3: Formation of FE equations in local
coordinates
W f   dTe NT fb dV   dTe NT f s dS  dTe (  NT fb dV )  dTe (  NT f s dS )
Ve
Se
Ve
Se
W f  dTe Fb  dTe Fs  dTe Fe
Fb 

N T f b dV
Fs 

NT f s dS
Se
Ve
fe  Fb  Fs
(Force vector)
Step 3: Formation of FE equations in local
coordinates


t2
t1
t2
t1
1T  1 T
( d e m ed e - d e k ed e  dTe Fe )dt  0
2
2
(Hamilton’s principle)
( dTe mede -  dTe k ede   dTe fe )dt  0
T
d
d
d
e

d   (
)  (dTe )
dt
dt
T
e


t2
t1
t2
t1
t2
t2
t2
T
T
T





 dt
d e m ed e dt  d e m e d e -  d e m ed e dt  -  dTe m ed
e
t1
t1
 - kd  F )dt  0 
dTe (-m ed
e
e
e
t1
  f
k ed e  m ed
e
e
FE Equation
Step 4: Coordinate transformation
  f
kde  m ed
e
e
(Local)
d e  TD e
y'
  F
K e De  M e D
e
e
(Global)
Global
coordinate
systems
x'
Local coordinate
systems
y'
y
where
x'
x
K e  T k e T , M e  T m e T , Fe  T f e
T
T
T
Step 5: Assembly of FE equations

Direct assembly method
– Adding up contributions made by elements
sharing the node
  F
KD  MD
KD  F
(Static)
  F
K e De  M e D
e
e
Step 6: Impose displacement constraints
No constraints  rigid body movement
(meaningless for static analysis)
 Remove rows and columns corresponding
to the degrees of freedom being constrained
 K is semi-positive definite

Step 7: Solve the FE equations

Solve the FE equation,
  F
KD  MD
for the displacement at the nodes, D

The strain and stress can be retrieved by
using e = LU and s = c e with the
interpolation, U=Nd
STATIC ANALYSIS

Solve KD=F for D
– Gauss elimination
– LU decomposition
– Etc.
EIGENVALUE ANALYSIS
  0
KD  MD
Assume
(Homogeneous equation, F = 0)
D   exp( it )
[K -  2M]  0
Let
 
2

[K - M]  0
det[K - M]  K - M  0
[ K - i M ]  i = 0
(Roots of equation are the
eigenvalues)
(Eigenvector)
EIGENVALUE ANALYSIS

Methods of solving eigenvalue equation
–
–
–
–
–
–
–
Jacobi’s method
Given’s method and Householder’s method
The bisection method (Sturm sequences)
Inverse iteration
QR method
Subspace iteration
Lanczos’ method
TRANSIENT ANALYSIS

Structure systems are very often subjected to
transient excitation.
 A transient excitation is a highly dynamic time
dependent force exerted on the structure, such as
earthquake, impact, and shocks.
 The discrete governing equation system usually
requires a different solver from that of eigenvalue
analysis.
 The widely used method is the so-called direct
integration method.
TRANSIENT ANALYSIS
(reading material)

The direct integration method is basically using
the finite difference method for time stepping.
 There are mainly two types of direct integration
method; one is implicit and the other is explicit.
 Implicit method (e.g. Newmark’s method) is more
efficient for relatively slow phenomena
 Explicit method (e.g. central differencing method)
is more efficient for very fast phenomena, such as
impact and explosion.
REMARKS

In FEM, the displacement field U is expressed by
displacements at nodes using shape functions N
defined over elements.
 The strain matrix B is the key in developing the
stiffness matrix.
 To develop FE equations for different types of
structure components, all that is needed to do is
define the shape function and then establish the
strain matrix B.
 The rest of the procedure is very much the same
for all types of elements.
Newmark’s method (Implicit)
Assume that

2  1

Dt t  Dt   t  Dt   t   -   Dt   Dt t 

 2

Dt t  Dt   t  1 -   Dt   Dt t 
 = 0.5
 = 0.25
Substitute into KD  CD  MD  F


2  1

K Dt   t  Dt   t   -   Dt   Dt t   

 2



Typically

C Dt   t  1 -   Dt   Dt t   MDt t  Ft t
Newmark’s method (Implicit)
K cm Dt t  Ftresidual
t
where
2
K cm  K  t   Ct  M 



21
 
Ftresidual

F
K
D


t
D


t






t
t t
t
t

 Dt  - C Dt   t 1 -   Dt 
2

 

-1 residual
Ft t
Therefore, Dt t  K cm
Newmark’s method (Implicit)
Start with D0 and D0
Obtain D0 using KD  CD  MD  F
-1 residual
cm t t
Obtain Dt using Dt t  K F
Obtain Dt and Dt using

2  1

Dt t  Dt   t  Dt   t   -   Dt   Dt t 

 2

Dt t  Dt   t  1 -   Dt   Dt t 
March
forward
in time
Central difference method (explicit)
MD  F - CD  KD   F - Fint  F residual
D  M -1F residual (Lumped mass – no need to solve matrix equation)
Dt t  2  t  Dt  Dt -t
Dt t  2  t  Dt  Dt -t
Dt 
1
 t 
2
 Dt t - 2Dt  Dt -t 
Dt -t  Dt -  t  Dt 
 t 
2
2
Dt
D0 and D0 are
prescribed and D0
can be obtained from
D,
D  M-1Fresidual
x
Central
difference
method (explicit)
x
x
Dt -t
x
x
Use Dt t / 2   t  Dt  Dt -t / 2 to
obtain Dt assuming Dt / 2  D0 .
Obtain Dt using D  M-1Fresidual
Obtain D-t using
2
t 

 Dt -  t  Dt 
Dt
2
t
-t -t/2 t0 t/2
Time marching in half the time step
Find average velocity D-t / 2 at time t =
-t/2 using D
t t / 2   t  Dt  Dt -t / 2
Find Dt / 2using the average acceleration at
time t = 0. D
  t  D  D
t t / 2
t
t -t / 2
Find Dt using the average velocity at time t =t/2
Dt t / 2   t  Dt  Dt -t / 2
t