Numerical Evaluation of

Dynamic Response
Dr Ignacio González

1
Outline
• Time stepping method’s
• Central Difference
• Newmark’s methods

• Duhamel’s integral

• (Repetition)

Ch. 5 ”Dynamic of Structures” A. Chopra

2
Reading in Chopra

• Duhamel’s integral (4.6-7)

• General about time stepping (5.1)
• Central Difference (5.3)
• Newmark’s Methods (5.4)
• STABILITY AND ERROR (5.5)

• Extra: Non-linear behavior (5.6-7)

3
In general, we want to solve the
”equation of motion”:

mu  cu  ku  p(t )

Given the initial conditions:
u0  u (0) u0  u (0)
Discretizing time and requiring the
equation of motion to be satisfied at a
finite number of points makes it a
computer-friendly problem.

4
 Continuous:
mu  cu  ku  p(t ) t tin  t  t fi

Discrete:

mui  cui  kui  pi i {0,1, 2 n}

We often only have time-discrete values of the load.

5
 Commonly the discretization points
equal are chosen equally spaced so
that:

ti 1  ti  t
ti
ti+1
Equilibrium equation are only
satisfied in the discretization
points. In between them
assumptions and simplifications
p

need to be made, introducing

t0 errors.

t

time

6
Every method should CONVERGE to the exact solution
when t  0.

There are two types of methods:

• EXPLICIT: use the equation of motion at t=ti plus
assumptions to find ui+1

• IMPLICIT: use the equation of motion at t=ti+1 plus

assumptions to find ui+1

7
Central Difference Method
(explicit)

Between any given three points ti-1, ti and ti+1 the movement is
assumed to be parabolic. This implies:

ui 1  ui 1 ui 1  2ui  ui 1
ui  ui 
2t t 2

mui  cui  kui  pi

8
 u i 1  u i 1 u i 1  2 u i  u i 1
u i  ui 
2t t2

m ui  c u i  k u i  p i

ui 1  2ui  ui 1 ui 1  ui 1
m c  kui  pi
t 2
2t

In the equation above only ui+1 solving it for that variable will
produce new values of u

9
To calculate u1 we need u0 and u-1, but what is u-1?
We know the position, speed and acceleration at time t=t0

u0 u0 initial conditions

p0  cu0  ku0
u0 
m
This leads to the 2 by 2 system below, where the unknowns
are u-1 and u1

u1  u1
u0 
t
u1  2u0  u1
u0 
t 2
10
Pseudo-code:

Initial calculations: For each step i:

p0  cu0  ku0
uo  pˆ i  pi  aui 1  bui
m
pˆ
t 2 ui 1  i
u1  u0  tu0  u0 kˆ
2
m c optional
kˆ  2  u u
t 2t ui  i 1 i 1
m c t
a 2 
t 2t u  2ui  ui 1
ui  i 1
2m t 2
bk 2
t

Check the book!!!

11
 Newmark’s Methods (implicit)
Between any given two points ti and ti+1 the acceleration is assumed to be
“something in between” the acceleration at ti and the acceleration at ti+1

ui 1  ui   a  t

ui 1  ui  1    ui     u11  t

ui 1  ui  ui t   a  t / 2
 2

ui 1  ui  ui t  1  2  ui   2  ui 1  t / 2
   2

12
 This leads to the 3 by 3 system above, with
unknowns ui+1, ύi+1,and ϋi+1

ui 1  ui  1    ui     ui 1  t

ui 1  ui  ui t  1  2  ui   2  ui 1  t / 2
2

mui 1  cui 1  kui 1  pi 1

Solving the system is left as homework for the

overambitious student

13
 Pseudo-code:
Initial calculations: For each step i:

pˆ i  pi  aui  bui

p0  cu0  ku0
u0  pˆ i
m ui 
 m kˆ
k k
ˆ c ui  ui   
t t 2 ui    t 1   ui
m  t   2 
a  c
t  ui ui ui
ui   
m    t t 2
2

b  t   1 c ui 1  ui  ui
2  2 
ui 1  ui  ui
ui 1  ui  ui

Check the book!!!

14
Explicit methods use the equilibrium at time ti and
“hope for the best”.

Implicit method use the equilibrium at time ti+1.

Therefore they can lead to wrong solutions but not

to “absurd” solutions.

This is generally not a problem since an accurate

solution needs a smaller step than a “not-absurd”
solution.

accurate inaccurate absurd

t
15
Harmonic driven damped
oscillator
• F=10 sin(/0.6 t) N
• m=0.2533 Kg
• k=10 N/m
• n=√ (k/m)=2 
• =5%;
• Ccr=2 mn c= 0.16 Ns/m
• u0=0;
• v0=0;

16
 Exact solution
t = 0.1 s Central Difference Method
t = 0.2 s (Explicit)
t = 0.3 s
t = 0.4 s
6
6

4 4

2 2

0 0

-2 -2

-4 -4

-6
0 0.5 1 1.5 2 2.5 3 -6
0 0.5 1 1.5 2 2.5 3

6
6

4
4

2
2

0
0

-2 -2

-4 -4

-6 -6
0 0.5 1 1.5 2 2.5 3 0 0.5 1 1.5 2 2.5 3

17
 Exact solution
t = 0.1 s Newmark’s =1/4 =1/2
t = 0.2 s (Implicit) “Average
t = 0.3 s Acceleration”
t = 0.4 s
6
5

0
0

-2

-4

-5 -6
0 0.5 1 1.5 2 2.5 3 3.5 4 4.5 5 0 0.5 1 1.5 2 2.5 3 3.5 4 4.5 5

0
0

-5 -5
0 0.5 1 1.5 2 2.5 3 3.5 4 4.5 5 0 0.5 1 1.5 2 2.5 3 3.5 4 4.5 5

18
 Exact solution
t = 0.1 s
Newmark’s =1/6 =1/2
t = 0.2 s (Implicit) “Linear
t = 0.3 s Acceleration”
t = 0.4 s
5 6

0 0

-2

-4

-5 -6
0 0.5 1 1.5 2 2.5 3 3.5 4 4.5 5 0 0.5 1 1.5 2 2.5 3 3.5 4 4.5 5
6
5
4

0
0
-2

-4

-6
-5
0 0.5 1 1.5 2 2.5 3 3.5 4 4.5 5 -8
0 0.5 1 1.5 2 2.5 3 3.5 4 4.5 5

19
 • Newmark’s linear acceleration
• Central difference
• Exact solution
Calculation time: 0.0007 s

3
displacement

-1

-2

-3

-4

-5
0 0.5 1 1.5 2 2.5 3 3.5 4 4.5 5

time

20
 Duhamel’s Integral

• An impulse can be seen as a spike-force of infinite

magnitude, zero duration and finite, non-zero integral over
time.
• The response to an impulse force can be exactly calculated.
• Any time-varying force can be seen as a sum of infinitely
many infinitely small impulses.
• Integrating the response to “each little impulse” we obtain
the exact solution for any time varying force.
• Errors are introduced if the integration is done numerically.

21
 d
Force is mass time acceleration: m u  p
dt
t2 t2
d
t1 p  dt  t1 m dt u  dt  m(u2  u1 )  m  u
It means an impulse is a change in the momentum,
a change in the velocity.
The response to an impulse is the response to a
“initial velocity”:

u0
u (t )  ent sin(d t )
d

h Impulse!
Where: u0 
m

22
Response to a single impulse h
applied at time ti:

h n (t ti )
u (t )  e sin(d (t  ti ))
md
Sum of N impulses aplied from ti
to ti+N df:

N
h n (t tik )
u (t )   e sin(d (t  ti  k ))
k  0 md

In the limit each “infinitely small” impulse is

h(t)=p(t) dt and the sum becomes and integral

p( )d n (t  )

t
u (t )   e sin(d (t   ))
t0
md

23
 Reordering:

t
1

n ( t  )
u (t )  e sin(d (t   )) p( )d
md t 0

This can’t be (easily) turned into a time-

stepping method. Why?

t 2 t t2


t1
f ( )d   f ( )d  f (t2 )t
t1

24
SDOF MDOF

Exact Mode Sup. Direct int.

Duhamel’s + dyn. factor

Explicit Implicit
Central diff. Newmark’s

25