SEISMIC / DYNAMIC

ANALYSIS USING FEM

Dhi
Dhiman
B
Basu

Acknowledgement
Text Book Dynamics
y
of Structures byy
Prof A K Chopra
y Text
e t Book
oo Dynamics
y a cs o
of Structures
St uctu es by
Prof J L Humar
y Text Book Finite
Finite Element Structural
Analysis by T Y Yang
y Most of the presentation materials are
taken from the course materials prepared
by Prof D C Rai,
Rai IIT Kanpur
y

REVIEW OF DYNAMICS AND

SEISMIC EXCITATION

Single Degree of Freedom (SDOF) System

One variable is required to describe the deformation
For example: Mass-Spring-Dashpot System

Newtons Law of Motion

Resultant force
= Mass X Acceleration

Resultant Force
Applied force, f (t )
Force offered by spring, kx
Force resisted by dashpot, cx

Equation of Motion

mx + cx + kx = f (t )

Second order, Non-homogeneous, Ordinary

Differential Equation
Initial Conditions:

Provided system parameters satisfy

certain condition

Multi-Degrees
Multi
Degrees of Freedom
5-Storied Building with
Columns axially inextensible
Beams are rigid
5 independent displacements are
required to specify the lateral
d f
deformation
i profile
fil
5 DOFs, one at each floor
5 (linearly) independent
deformation shape
Mode-1
5
4
3
2
1
0

Mode-2

Mode-3

5
5
4
4
3
3
2
2
1
1
0
0
0 0.25 0.5 0.75 1 -1.5 -1 -0.5 0 0.5 1 -1.5 -1 -0.5 0 0.5 1

Mode-4

Mode-5

5
4
3
2
1
0
-2 -1.5 -1 -0.5 0 0.5 1

5
4
3
2
1
0
-4 -3 -2 -1 0 1 2 3

Any instantaneous deformation profile is a linear combination of these

shapes

Each storey can be modeled

with a mass- spring-dashpot
assembly
Floor mass contributes to the
mass
Columns contribute to the
spring stiffness
Dashpot accounts for the
inherent damping
NOTE the eq
equivalence
i alence of 2nd and 3rd figures
fi res Springs
Sprin s
and dashpots are acting laterally

Dynamic Equilibrium
Dynamic Equilibrium of MDOF
System
[ M ]{u} + [C ]{u } + [ K ]{u } = {P}
[ M ] = Mass Matrix (Diagonal here)
[C ] = Damping Matrix (Tri-diagonal here)
[ K ] = Stiffness Matrix (Tri-diagonal here)
{P} = Excitation Vector

[M] (and [C]) has the similar

definition as that of [K]

Force Equilibrium
Applied Force=Resistive Force

Applied Force = Pi (t )

Arbitrarily time varying force

Under Seismic Excitation

Spatially Uniform
Ground Motion

[ M ]{u} + [C ]{u } + [ K ]{u } = {P}

{P} = [ M ]{A}ug

Characterization of Ground Motion

Accelerration (g)

0.2
0.1

EW

0
-0.1
-0.2
0

10

Accceleration (g)

02
0.2

20
Time (sec)

30

40

0.1

NS

0
01
-0.1
-0.2
0

10

20
Time (sec)

30

40

Acceleration (g)
A

0.2
0.1

Vertical

0
-0.1
-0.2
0

10

20
Time (sec)

30

40

mx + cx + kx = f (t ) = mug (t )

Spectral acce
S
eleration (g)

Response Spectrum
06
0.6

x + 2 n x + n2 x = ug (t )

0.4

EW

S = MaxR(Tn , )

0.2

Spectral acce
S
eleration (g)

2
4
Period (sec)

06
0.6

0.4

NS
Vertical

0.2

0
0

2
4
Period (sec)

REVIEW OF ANALYTICAL MECHANICS

Generalized Coordinate and

Constraint Equations
yi = fi (q1 , q2 ,...)
yi
qjs
j=1 q j
N

yi =

Virtual Displacement
Work done = ZERO
M

i=1

i=1

yi
qj

q
j=1
j
N

Fi yi = Fi

N
M
yi

q j = Q j q j
= Fi
i=1 q j
j=1
j=1
N

Generalized Force
M

Q j = Fi
i=1

yi
= zero for Equilibrium
q j


L
y1 = q1 q2 j

L

y2 = q1 a q2 j

Relation between Generalized and

Constraint Coordinate

y3 = (q1 ) j

L
y4 = q1 + q2 j

L
(k1 k2 ) P mg
k1 + k2

2
q1
=

L k k

q2
P a 2
(
)

1
2
2

Equilibrium
E
ilib i
Equation
E
i
obtained by equating
generalized forces to zero

Constraint Equations
g1 ( y1 , y2 ,..., t ) = 0
g 2 ( y1 , y2 ,..., t ) = 0

Holonomic constraints

Example
p of Constraint in Building
g Analysis
y

3 Dofs Per Floor Level

Lagranges
Lagrange s Equations
In Generalized Coordinate

T = T (q1 , q2 ,...., q1 , q 2 ,....)

Kinetic Energy

V = V (q1 , q2 ,....)

P tential Ener
Potential
Energy

Lagranges Equation
d T T V

+
Qj = 0

dt q j q j q j
Qj =

Those forces that cannot be derived from a scalar

function, e.g., damping force, externally applied forces
etc.

STATIC ANALYSIS TO DYNAMIC

ANALYSIS

Vibrating Flexural Beam Element

Potential Energy V = Strain Energy U =
EI
2

= mass density
d i
A= area
Lagrange Equation

2 v
x 2 dx
L

Kinetic
i i Energy T =

A
2

v ( x ) dx

d T T V

+
Qj = 0



q
q
q

dt
j
j
j

L
L

2 v 2 v

d T V
d A

dx
Y1 = +
v ( x )dx + EI 2
=
2 v ( x )

v1 x 2

dt v1 v1 dt 2 0
x
v1

= (m11v1 +m121 +m13v 2 +m14 2 ) + ( K11v 1 + K12 1 + K13v 2 + K14 2 )

L

K ij = EI fi '' ( x) f j'' ( x)dx

0
L

mij = A
0

f i ( x ) f j ( x)dx

 Y1
156
22 L

M 1 AL 22 L
4 L2

=
Y2 420 54
13L

13L 3L2
M 2

12
6
12
6
2

2
L
L
L
L

v
54
13L v1
6
6

1

4
2

13L 3L 1 EI L
L
1
+

156 22 L v2 L 12
6 12
6 v2


2
2
2
L
L
L
L 2
22 L 4 L 2

2
4

L
L

Consistent Mass Matrix

Damping Terms can be added separately
Element Equilibrium Equations can be assembled to form
Global structural equilibrium
q

USE of FEM in DYNAMIC ANALYSIS

Overview of FEM for Static Loads

Divide the continuum into a finite number of subregions

g
((elements)) of simple
p g
geometry
y
Select key points on the elements to serve as nodes,
where conditions of equilibrium and compatibility are to
be enforced.
Satisfy strain-displacement
strain displacement and stress-strain
stress strain
dy
relationships within a typical element
Determine stiffnesses and equivalent nodal loads for a
typical element using work or energy principles
Develop equilibrium equations for the node of
discretized continuum in terms of the element
K {u}
contributions
Solve these equilibrium equations for the nodal
displacements
Calculate stresses at selected points within the elements
Determine support reactions at restrained nodes if
d i d
desired

dx

dz

{u}

Overview of FEM for Dynamic

Loads

Divide the continuum into a finite number of sub-regions of

simple geometry
Select key points on the elements to serve as nodes, where
conditions of equilibrium and compatibility are to be enforced.
Satisfy strain
strain-displacement
displacement and stress
stress-strain
strain relationships within
a typical element
Determine stiffness, mass, damping and equivalent nodal
dynamic loads for a typical element using work or energy
d
dy
principles
Develop dynamic equilibrium equations for the node of
discretized continuum in terms of the element contribution
S l these
Solve
h
equilibrium
lb
equations ffor the
h nodal
d l displacements
d l
Calculate stresses at selected points within the elements
Determine support reactions at restrained nodes if desired

dV {u}
cdV {u}

dx

dz

{u}

FEM for Dynamic Loads

Displacement formulation
Displacement field is function of time

{u } = {u
f

( x, y , z , t )

Nodal displacements are functions of time

{u} = {u (t )}
The spatial interpolation is unchanged

u f ( x, y, z , t ) = [ N ( x, y, z ) ]{u (t )}

N Shape function

dy

dx

dz

Equations of Motion
dy

dx

dz

Infinitesimal element in 3D
The
Th spatial
ti l iinterpolation
t
l ti

dV {u}
cdV {u}

{u}

For general volumetric loads

f e.n.l = N T ( cu u) dV





V
f e.n.l = cN T NdVu N T NdVu
   V  

V
f e.n.l = cu mu
  


{u} = N ( x, y, z ) {u ( t )}
{u} = N ( x, y, z ) {u ( t )} Apply DAlemberts principle
ku = f e.n.l + p (t )
{u} = N ( x, y, z ) {u( t )}
 
mu + cu + ku = p (t )
   

Force and Base Excitations

Equation of motion for system is

mu + cu + ku = p (t )
    

Special case is seismic load in 3-Dimension

mu + cu + ku = mx uxg m y uyg mz uzg

   

dy

dx

dz

p(t)

Resultant force vector

p (t ) = mx uxg m y uyg mz uzg

dy

dx

dz

uigg
uig is the ground acceleration in direction i and m is a column matrix which
i

represents the sum of all columns in the mass matrix m associated with
displacements in the i-direction
This definition is valid only if the vector is defined as the displacement relative to
the displacement at the base of the structure

Overview of Dynamic Problems and

Analysis Types
Problems

Analysis

Rigid body vs elastic

Real eigen value
motion
analysis
Eigen value
Free vibration vs forced Complex eigen value
vibration
analysis
Loading
g type
yp
Frequency
q
y response
p
Response
o Harmonic
Transient response
o Periodic
Random response
o
o
o

Transient
T
i
Stationary random
random

FE Modeling for Dynamic Analysis

Continuous Model
Infinite degree of freedom
Equation of motion is PDE
Finite Element Model
Continuous mass is discretized
into nodal masses
Finite degree of freedom
Assumed displacement field
Equation
E ti off motion
ti iis ODE
Discretization of mass can be done
usingg two formulations
Consistent Mass
Lumped Mass

PDE

x
u(x,t)
( t)
x1

ODE

u(x1,t)

ODE

x1

x2

u(x1,t)

u(x2,t)
t)

x3

u(x3,t)

Consistent Mass Matrix

General approach for element mass matrix
m = N T NdV
 V  

Alternativelyy

m = h 1 g T gdVh 1

 
V

u f = gc


u = hc
   1

N = gh


where g = geometric matrix


c
= undermined
d
i d constants
t t
c = h 1u

f
 f   1
u=
displacement field

u = gh u
u= nodal displacement

 

N= shape function


For different types of elements, geometric matrix is different and

hence equivalent mass matrix is also different

For 1D flexure beam element, considering the

translational inertia only
m = h 1 g T gdVh 1

 
V
g = 1 x

1

x
T
g g= 2
  x
3
x

x3

x2
x

x2
x3
x4

x3
x5
x5

x4
x6

x 6
3

1 L2
L
2
T
g
gdx
=
0   1 3
L
3
1
L4
4

Consistent Element Mass Matrix

156 22 L

2
AL 22 L 4 L
m=
13L
420 54


2
13L 3L

1 2
L
2
1 3
L
3
1 4
L
4
1 5
L
5

1 3
L
3
1 4
L
4
1 5
L
5
1 6
L
6

54
13L
13L 3L2
156 22 L

22 L 4 L2

1 4
L
4

1 5
L
5 1
h
1 6 
L
6
1 7
L
7

L3
0

L3
1 0
= 3
L 3L 2 L2

L
2

Beam
Element

0 0
3L L2

2 L
0

For 1D flexure beam element,, considering

g the rotational
inertia only

mr = h 1 g, xT g, x dVh 1

 
V
0 0
0 1
g, xT g, x =
0 2 x
 

2
0 3x

0
2x
4 x2
6 x3

where

g, x = 0 1 2 x 3 x 2


0 0
0 L

T
2
=
g
g
dx
x
x
,
,
0  
0 L

0 L3

0
3 x 2
6 x3

9x4

36
Rotational Mass

Matrix for
I 3L
r =
Beam Element m
30 L 36


L2
4 3
L
3
3 4
L
2

0
L3
3 4
L
2
9 5
L
5

36
3 L

3L
L2
3L 36 3L

2
2

3
L

3
L
4
L

3L
4 L2

Beam
Element

Assemblage of Property Matrices

Two beam elements with three 1
unrestrained nodes

Assembly of stiffness matrix

12
6L

E I 12
[K ] = 3 6 L
L
0

6L
4 L2
6 L
2 L2

12
6 L
12 + 12
6 L + 6 L

6L
2 L2
6 L + 6 L
4 L 2 + 4 L2

0
0
12
6 L

0
0

12
6L

6 L
2 L2

12
6 L

6 L

4 L2
0
0
6L
2 L2

2
L

1
2
3
4
5
6

22 L

54

13 L

4L

13 L

3L

13 L

156 + 156

22 L + 22 L

22 L

3L
0

22 L + 22 L
54

4L + 4L
13 L

13 L
156

13 L

3 L2

22 L

0
13 L

3 L2
22 L

4 L2
0

1
2
3
4
5
6

3
L

Similarly, Assembly of Mass (consistent) matrix

156
22 L

AL 54
[M ] =

420 13 L
0

Applying Boundary Conditions

Exclude rigid body modes (typical for civil engineering problems)
For example, in two element cantilever beam, Node 3 is restrained node
so DOFs 5 & 6 can be eliminated

Mass matrix for constrained beam

156 22 L 54 13L

2
2
AL 22 L 4 L 13L 3L
m=
13L 312
0
 420 54

2
2
L
L
L

13
3
0
8

1
2
3
4

L
L1

Stiffness matrix for constrained beam

6 L 12 6 L
12

2
2
EI 6 L 4 L 6 L 2 L
k= 3
0L
 L 12 6 L 24

2
2
6
L
2
L
0
8
L

1
2
3
4

Reduction of Size for Dynamic

Analysis

Stiffness (static) reduction

mRR

m
TR


mRT uR k RR
 + 
mTT uT kTR



k RR uR + k RT uT = 0

1
uR = k RR k RT uT
 

k RT uR FR = 0
 =
kTT uT FT

1
3
2

Hence, from the second equation

TT T

k u = FT
 
where

k TT
= k RR kTR k RR1 k RT = 0


  

Mass (dynamic) reduction

1
k RR

k
RT
m = T mTT TT =  
 

 IT

Applying stiffness reduction

TT

kRR
k
TR


T
T

4L2 2L2 6L 6L

kRT EI 2L2 8L2 6L

0
=

kTT L3 6L 6L 12 12


6
0
12
24

EI 1.7143

k TT
= 3

L 4.2857

4
1

4.2857
13.714

1
3

Similarly, applying mass reduction

mRR

m
TR


4L2 3L2 22L 13L 2

mRT AL 3L2 8L2 13L 0 4

 =
mTT 420 22L 13L 156 54 1


L
13
0
54
312

1
0.27464
027464
0.17566
017566

mTT
= AL
3
0.17566
0.89096


1
3

Vibration of a Cantilever
Free vibration of a uniform cantilever beam
L

4u
2u
in t
L=1000mm, A=(20x5)mm2
u
x
,
t
=

x
e
EI 4 + AL 2 = 0 ( ) n ( )
E = 2.1x105 N/mm2
x
t
n ( x ) = characteristic function (nth mode)
n4 = AL n2 / EI
d 4n ( x )
4

n ( x ) = 0
n
4
n = ( n l 2 ) EI / m 4
dx
General solution is

= natural frequency of nth mode

n ( x ) = c1 sinh x + c2 cosh x + c3 sinh x + c4 cosh x

Applying BCs and solving the characteristic equation natural frequencies of
beam can be obtained

1
EI
3.156 4
mL

2
1

EI
22.03 4
mL

3
1

EI
61.70 4
mL

Equation of motion for undamped system is

mu + Ku = 0
 
u (t ) = cos t

k 2 m = 0

 

Using k TT
and mTT
the eigen value problem is setup


k11 i 2 m11 k12 i 2 m12 1,i 0

=
2
2
(i = 1
1, 2)
k21 i m21 k22 i m22 2,2 i 0

Solving above equation and substituting L = L1/2

1 = 3.5222

1
L12

1
2 = 22.28
22 28 2
L1

EI
A
EI
A

= [ 1

1.0000 1.0000
2 ] =

0.3395 0.6991

Lumped Mass Matrix

Simpler model for inertia properties for

obtaining diagonal mass matrices

Various approaches

(1)
m1

Tributary masses are lumped at the nodes of

discretized continuum
Lumped mass matrix

m1
0

...
m=
0
...

(2)
m2
L1

m( 2)/2
m( 1)/2

...

...

m2
...
0
...
0

... 0 ...
... ... ...
... mi ...
... ... ...
... 0 ...

0
0
...

0
...

mnn

1
mi = ( m(i ) + m(i +1) )
2
where nn = number of nodes
mi = tributary masses lumped at node i

Using the same k TT and Lumped Mass Matrix (LMM)mLTT



different set of and are obtained

m1
L

m2

L1

k11 i 2 m11

2
k

i m21
21

mTT


0.5 0
= AL

0
1.0

1
3

k12 i 2 m12 1,i 0

= (i = 1, 2)
2
k22 i m22 2,i 0

Solving above equation and substituting L = L1/2

1
1 = 3.156 2
L1
1
2 = 16.2584 2
L1

EI
A
EI
A

= [ 1

1.0000 1.0000
2 ] =

0.3273 1.5272

Comparison of Eigen Soultions

Consistent Mass Matrix
1 = 3.5222
2 = 22.28

1
L12

1
L12

Lumped Mass Matrix

EI
A

1 = 3.156

EI
A

2 = 16.2584

1.0000 1.0000
= [ 1 2 ] =
0 3395 00.6991
6991
0.3395

EI
A

1
L12
1
L12

EI
A

1.0000 1.0000
= [ 1 2 ] =
0 3273 11.5272
5272
0.3273
(0.339)CM
(0.327)LM

Distributed mass model

1
EI
3 156 4
3.156
mL

(1.0)CM,LM

2
1

EI
22 03 4
22.03
mL

(-1.53)LM

(-0.699)LM

(1.0)CM,LM

Effect of Discretization
Lumped mass matrix (LPM)

L=1000mm,
A=(20x5)mm2
E = 2.1x10^5 N/mm2

More number of elements

better solution

Consistent mass matrix (CMM)

LPM -lower bound solution,

CMM - upper bound solution
For same number of element
CMM shows better results

Structural Damping Matrix

Proportionally
p
y Damped
p System
y
Equation of motion can be dynamically uncoupled by the eigenvectors of the
undamped system

mu + cP u + ku = 0
    
Rayleigh Damping

cP = m + k




2i j

i + j

c = m+k

c = k

2
=
i + j

= damping ratio
i j = natural ffrequency
q
y off ith and jth mode

c = m
r

Transform to normal coordinates

[CP ] = T cP = 2 rr M r


Response History Analysis

Governing Equation

[ M ]{u} + [C ]{u } + [ K ]{u } = {P}

{P} = [ M ]{A}ug

Analysis of Building using SAP2000

Ground Motion
0.2

Acc
celeration (g)

0.15
0.1
0.05
0
-0.05

10

15

20

25

30

35

-0.1
-0.15
-0.2

Time (sec)

Ground Motion Along Y-direction

40

Base Shear

In KN

Roof Displacement at One Corner

In mm

THANK YOU