Structural Dynamics

Structural dynamics 1.1

SDF - INTRODUCTION A single degree of freedom system consisting of a mass m and a spring with stiffness k is considered.

u(t ) k m

k (N/m)

p(t )

The spring is undeformed for u = 0

D’Alembert’s principle The system is supposed in dynamic equilibrium. The principles of statics are applied by introducing a fictitious inertia force, a force equal to the product of mass time its acceleration and acting in a direction opposite to the acceleration.

Free-body diagram

The equation of motion can be derived in two ways.

u&&(t ) ku(t )

p(t )

Newton’s second law mu&&(t )

u&&(t )

u(t ) > 0 k m

p(t )

m u&&(t ) = p(t ) − k u(t )

ku(t )

fictitious inertia force p(t )

Dynamic equilibrium →

m u&& + k u = p(t )

p(t ) − k u(t ) − m u&&(t ) = 0


Stiffness k Combination of two springs – case 2

The stiffness k has the same definition as in the displacement method: The stiffness k is the external force that is needed to keep the system in equilibrium when a unit displacement u = 1 is applied.

spring

m

k2

This case is more complicated.

m

u =1

u =1

k.1

k1

F1 = k1 d1

d1, F1 d 2 , F2 m

k

k

F2 = k 2 d 2

external force

statics : F1 = F2 = k Combination of two springs – case 1

k2 k1

The stiffness for this system is trivial:

d1 + d2 = 1 →

k k + =1 k1 k 2

m

k = k1 + k 2

1 1 1 = + k k1 k 2


Some structures can be idealised as SDF

Example

3

2

5

1

6

The mass of this SDF system is m, the mass of the roof.

4

In statics, this frame has 6 active degrees of freedom.

The stiffness is determined in the classical way: u =1 2

3

k

1

k=

By neglecting the axial deformations, 3 d.o.f. disappear.

1

Only one d.o.f. is left if the frame is consisting of an heavy roof supported by light columns.

m

massless

rigid beam

12EI 3

1

6EI 2

L

L

6EI

12EI

2

L

L

3

24EI L3


Rotations

A SDF system can also have a rotational movement.

L

d

Jθ&& is the moment created by the fictitious inertia forces.

massless

G

M(t)

m

θ

Example: Calculation of J for a bar

R

dm x θ&&

in statics :

M =Cθ

π d4 G C= 32 L

inertia force fI m

ρ=

θ

x

m L

L

in dynamics :

J θ&& + C θ = M (t ) L

Jθ&& = ∑ fI ⋅ x = ∫ x 2 d m ⋅ θ&& 0

J is the moment of inertia . for a circular disk :

J=

mR 2

2

L3 mL2 = J = ∫ x dm = ρ ∫ x d x = ρ 3 3 0 0 L

2

L

2


SDF - FREE VIBRATION u(t ) = uo cos(ωnt ) +

UNDAMPED FREE VIBRATION The structure is disturbed from its static equilibrium and then vibrates without any applied forces. The equation of motion is:

The solution is:

ωn = k m

u(t ) which can be written as

k

ωn

sin(ωnt )

u(t ) = C sin(ωnt + θ )

m

C = uo2 + (u&o ωn )2 m u&& + k u = 0

u(t ) = A cos(ωnt ) + B sin(ωnt ) (rad/s)

u&o

cosθ =

u&o ωn C

sinθ =

u

u&o

C

uo

t

natural circular frequency Tn

A and B are determined by the initial conditions natural period

ut = 0 = uo → uo = A u&t = 0 = u&o → u&o = B ωn

uo C

Tn =

2π

ωn

(s)

natural frequency

fn =

1 ωn = (Hz) Tn 2π


Energy in undamped free vibration

At each instant of time, the total energy E is made of two parts, the kinetic energy Ek and the strain energy Es.

u E = max s E =0 k

t

Es (t ) =

1 k u(t )2 2

Ek (t ) =

1 m u& (t )2 2

E = max s E =0 k

E =0 s E = max k

E (t ) = Ek (t ) + Es (t ) 1 1 2 2 = m [C ωn cos(ωnt + θ )] + k [C sin(ωnt + θ )] 2 2 1 1 = m C2 ω n2 cos2 (ω nt + θ ) + k C2 sin2 (ωnt + θ ) 2 2 1 = k C2 2

( k = mωn2 )

E(t) is constant, which implies conservation of energy.

Remark : the conservation of energy can be used to derive the differential equation.

1 1 E (t ) = m u& 2 + k u 2 2 2

conservation of energy → dE = mu&u&& + kuu& = 0 dt

→

dE =0 dt

mu&& + ku = 0


VISCOUSLY DAMPED FREE VIBRATION A and B are determined by the initial conditions. Friction in the structure is idealized by a linear viscous damper which develops a c force proportional to the k velocity fD = − c u& (t ) The equation of motion is:

u(t ) u& (t ) = −ξωn e −ξωnt [A cos(ωDt ) + B sin(ω Dt )] m

m u&& + c u& + k u = 0

if c < cr = 2 km (critical damping) the solution is u(t ) = e −ξωnt [A cos(ωDt ) + B sin(ωDt )]

damping ratio

ξ=

c c = cr 2 km

damped pulsation

ωD = ωn 1 − ξ 2

+ e −ξωnt [− A ωD sin(ω Dt ) + B ω D cos(ωDt )]

ut =0 = uo

→

uo = A

u&t =0 = u&o

→

u&o = −ξωn A + B ωD

  u& + ξωnuo u(t ) = e −ξωnt uo cos(ωDt ) + o sin(ω Dt ) ωD   The solution can also be written as

u(t ) = C e −ξωnt sin(ωDt + θ )  u& + ξωnuo  C = uo2 +  o  ω   D

2

cos θ =

sinθ =

uo C

u&o + ξωnuo C ωD


u

C e −ξωnt

p periods between two maximal points un and un + p

un

tn + p = tn + pTD

un + p t n+p

tn

TD =

2π

t

C e − ξωntn sin( ωDtn + θ ) un = un + p C e − ξωn (t n + pTD ) sin( ωD (t n + pTD ) + θ ) = eξωn pTD

ωD

un 2π = ξ ωn p TD = ξ ωn p un + p ωn 1 − ξ 2

Decay of motion

ln

A free vibration test can be used to determine experimentally the natural frequency and the damping of a structure.

ξ < 0 .1 →

1−ξ2 ≈1 → ξ =

1 u ln n 2 πp u n + p


COULOMB-DAMPED FREE VIBRATION Coulomb damping results from friction against sliding of two dry surfaces.

u(t ) k m

The constants A1, B1, A2, B2 depend on the initial conditions of each successive half-cycle motion.

dry friction

The friction force is F = µ N where µ denotes the coefficients of static and kinetic friction, taken to be equal, and N the normal force between the sliding surfaces.

The plot of the solution is u

F is assumed to be independent of the velocity of the motion and its direction opposes motion. The equations of motion from left to right are m u&& + k u = −F

u(t ) = A1 cos(ωnt ) + B1 sin(ωnt ) − F / k The equations of motion from right to left are m u&& + k u = F

u(t ) = A 2 cos(ωnt ) + B2 sin(ωnt ) + F / k

t


DAMPING MULTIPLE DEGREES OF FREEDOM Different damping models can be used, a viscous proportional damping is the most used approach.

There are two reasons for that: • The mathematical equation which describes the motion is easy.

SDF

m u&& + c u& + k u = p(t )

MDF

[ m] {u&&} + [c] {u& } + [ k ] { u} = {p(t )} [ k]

stiffnesss matrix

[ c]

damping matrix

[ m]

mass matrix

• This model gives results which are often in very good agreement with experiments.

A consequence is that the damping coefficient ξ can only be determined by experiments.

The mass and stiffness matrices are obtained by (finite element) discretisation of the structure. The damping matrix cannot be obtained by discretisation, a different approach must be used.


SDF - HARMONIC EXCITATIONS Example 2 A harmonic load is applied to the structure.

y

u(t ) k

u po sin(ωt )

m

Without damping, the equation of motion is

y = u + ug

m y&& = −k u

k /2

k /2

ug (t ) = ugo sin(ωt )

m u&& + k u = po sin(ωt ) Example 1 The system consisting of the mass m and the excentric mass m is considered for writing Newton’s equation.

m

If u is to be studied (e.g. earthquake)

m ( u&& + u&&g ) = − k u

m

v = r sin(ωt )

→

m

m u&& + k u = m ugo ω 2 sin(ωt )

u

k /2 k /2

(m − m ) u&& + m ( v&& + u&& ) = − k u m u&& + k u = m r ω 2 sin(ωt )

If y is to be studied (e.g. floor isolation)

m y&& = − k ( y − ug ) →

m y&& + k y = k ugo sin(ωt )


WITHOUT DAMPING

C and θ are determined by the initial conditions

u(t )

po sin(ωt )

k

m u&& + k u = po sin(ωt )

m

The solution u(t) of the differential equation is the sum of two parts uh (t) and up (t).

u(t ) = uh (t ) + u p (t )

u(t) is a summation of two sinus and is not defined for ω = ωn. In reality, the damping implies that uh(t) disappears after some time. Then the solution ( steady state response ) is

u(t ) = u p (t ) =

homogeneous solution uh (t ) = C sin(ωnt + θ ) particular solution

u&&p (t ) = − Aω 2 sin(ωt )

u(t ) = C sin(ω nt + θ ) +

po k − mω2

sin(ωt )

po / k 1 − (ω ωn )2

• After some time, the structure vibrates with the same frequency as the applied force. • The amplitudes of the vibration are infinite when ω = ωn. uh(t) is the transient response

po k 1 − (ω ω n )

=

1 − (ω ωn )

2

conclusions

u p (t ) = A sin(ωt )

→ − m A ω 2 + kA = po → A =

po k

2

sin(ωt )

up(t) is the steady state response


WITH VISCOUS DAMPING after calculations (see the book), it is obtained

u(t ) c k

po sin(ωt ) m

po k

u p (t ) =

[1 − (ω ω ) ] n

m u&& + c u& + k u = po sin(ωt ) The homogeneous solution uh(t) (transient response) disappears after some time.

uh (t ) = C e−ξωnt sin(ω Dt + θ ) The particular solution up(t) (steady state response) is of the form:

u p (t ) = A sin(ωt − φ )

tanφ =

2 2

+ [2ξ (ω ωn )]

2ξ (ω ωn ) 1 − (ω ω n )2

sin(ωt − φ ) 2

0 < φ < 180o

Remarks • The total response is u(t) = uh(t) + up(t). But after some time uh(t) disappears and u(t) = up(t) (steady state response). • After some time, the structure vibrates with the same frequency as the applied force.


This numerical example shows that the transient response uh(t) disappears after some time and that only the steady state response up(t) is then left.

1 0.8 0.6

u 0.4 0.2

u(t ) = uh (t ) + u p (t )

0 -0.2 -0.4

u(t ) = C e −ξωnt sin(ω Dt + θ ) +

-0.6

po k

[1 − (ω ω ) ] n

2 2

+ [2ξ (ω ωn )]

sin(ωt − φ ) 2

-0.8

0

5

10

15

20

t total response u(t) transient response uh(t)

C and θ are determined by the initial conditions. For this case, the steady state is obtained after about 20 seconds.

25


Dynamic factor After some while, the structure vibrates with the same frequency as the applied force. It is the steady state response up(t).

Rd can be plotted as function of the ratio ω /ωn for different values of the damping coefficient ξ. 5 4.5

The amplitude of these vibrations are now studied.

ξ = 0.1

4 3.5

u p (t ) =

po k

[1 − (ω ω ) ] n

2 2

+ [2ξ (ω ωn )]

sin(ωt − φ ) 2

Rd

3

ξ = 0.2

2.5 2

ξ = 0.4

1.5

The static deformation due to a static load po is

1

p (ust )o = o k

0.5 0

The amplitude of the vibration is equal to the product of the static deformation times a dimensionless dynamic factor Rd.

0

0.5

1

1.5

2

2.5

3

(ω ω n )

ω ω n < 0.25 → Rd ≈ 1 " quasi static" response m u&& + c u& + k u = po sin(ωt ) → u = po k sin(ωt )

Rd (ξ ,ω ωn ) =

[1 − (ω ω ) ] n

1

2 2

+ [2ξ (ω ωn )]

2

ω → ωn the amplitudes of vibrations become large : Resonance


Resonance

Band-width method ( Experimental method to determine ξ )

Resonance is reached for

ω = ωr = ωn 1 − 2ξ 2 For this value of ω , the dynamic factor is

Rd max =

1

2ξ 1 − ξ 2

if ξ < 0.1 then Rd max ≈

1 and 2ξ

The structure is excited by a harmonic load. The frequency of the load is increased step by step. At each step, the amplitudes of vibrations of the steady state response are measured. This implies that at each step, some time must be waited so that the transient response disappears. The curve Rd as function of ω is then obtained experimentally. 3

Rd max

Rd 2

ωn = k m ≈ ω D = ωn 1 − ξ 2 ≈ ωr = ωn 1 − 2ξ 2

1

Example

0

0

ξ = 2% → Rd = 25 The deformations are 25 times the static ones.

ξ < 0.1 → ξ ≈

Rd max 2

ω1 ω r ω 2

ω2 − ω1 f2 − f1 = ω 2 + ω1 f2 + f1

ω

2

f =

ω 2π


Force transmission and vibration isolation A harmonic load is applied to a structure. This structure is connected to the ground through a support modelled by a spring k and a damper c.

p(t ) = po sin(ωt ) u Steady state response

m k

c

u (t ) =

fT (t ) = poRd [sin(ωt − φ ) + 2ξ (ω ωn ) cos(ωt − φ )] = poRd 1 + [2ξ (ω ωn )] sin(ωt − φ + α ) 2

The transmissibility TR is defined as the ratio between the amplitude of the transmitted force fT and the amplitude of force applied to the structure.

po Rd sin(ωt − φ ) k f TR = T max = po

fT

1 + [2ξ (ω ωn )]

[1 − (ω ω ) ]

2 2

n

2

+ [2ξ (ω ωn )]

2

The force transmitted to the ground is TR is dimensionless.

fT (t ) = ku(t ) + cu& (t ) = poRd sin(ωt − φ ) +

poc Rdω cos(ωt − φ ) k

The objective is to choose the support (k, c) such that TR is as small as possible.


The mass m of the structure and the frequency of the load ω are imposed. The problem is to choose a support with k and ξ such that the transmissibility TR is minimal. TR as function of the ratio ω /ωn is plotted for different values of ξ. ωn = k m

3

ξ =0 2.5

ξ = 0.2 ξ = 0.25

TR The transmitted force is less than the applied one if the natural frequency ωn is such that ω /ωn > 1.4 . A low TR is obtained for low values of ωn and ξ. However, a very low ωn implies a low k and therefore a too large static displacement po / k. Besides, a very low ξ implies high displacement amplitude while passing through resonance which may occur before the load reaches the circular frequency ω. Natural rubber is a good compromise and is often used for the isolation of vibrations.

2

1.5

ξ = 0.5

1

0.5

0 0

0.5

1

1.5

ω ωn

2

2.5

3


SDF - ARBITRARY EXCITATIONS 0 ≤ t ≤ to

m u&& + k u = po

This chapter studies the response of a SDF system to pulse, impulse and periodic excitations.

u(t ) = uh (t ) + u p (t )

RECTANGULAR PULSE FORCE uh (t ) = A cos(ωnt ) + B sin(ωnt ) p(t )

u p (t ) =

u(t ) k m

p(t )

po

to initial conditions:

t

uo = u&o = 0

An undamped SDF system is loaded by a rectangular pulse force. For this case, the differential equation can be solved analytically.

po k

u(t ) = A cos(ωnt ) + B sin(ωnt ) +

 uo = 0    u&o = 0

→

u(t ) =

po k

po  A + =0  k    Bωn = 0

po [1 − cos(ωnt )] k


t ≥ to

m u&& + k u = 0

u(t ) = u(to ) cos(ωnt1 ) +

t1 = t − to

free vibration

u& (to )

ωn

sin(ωnt1 )

po [cos(ωnt − ωnto ) − cos(ωnto ) cos(ωnt − ωnto ) k

+ sin(ωnto ) sin(ωnt − ωnto )]

u(to ) =

po [1 − cos(ωnto )] k

by using

u& (to ) =

poωn sin(ωnto ) k

it is obtained

p u(t ) = o [1 − cos(ωnto )] cos(ωnt − ωnto ) k +

u(t ) =

po sin(ωnto ) sin(ωnt − ωnto ) k

u(t ) =

cosa cosb − sina sin b = cos(a + b)

po [cos(ωnt − ωnto ) − cos(ωnt )] k


IMPULSE LOAD p(t )

uo = u&o = 0

1 d u& (t ) = p(t ) d t m

u(t )

→

k m

p(t )

to

t

An undamped SDF system is loaded by a short impulse force. An approximative response is to be calculated. Newton’s law for t ≤ to

d u& (t ) m + k u(t ) = p(t ) dt

Hypothesis: to is so small that the displacement u is still zero at to. (uo = 0).

 to is very small → u(to ) ≈ 0  u = 0  o Then ku(t) can be neglected in Newton’s equation

u& (to )

→ u& (to ) − u&o

1 to = ∫ p(t ) d t m0 to

I = m

I = ∫ p(t ) d t 0

At to, the impulse load is assumed to have produced an initial velocity, but no displacement. Newton’s law for t ≥ to

(free vibration)

u(t ) = u(to ) cos[ωn (t − to )] +

u& (to )

ωn

sin[ωn (t − to )]

→ the maximum displaceme nt is umax =

I mω n

This result is only valid if to is small enough so that u(to) ≈ 0. In practice it means to < Tn /10 In such case, the maximal deformation does not depend on the form of the impulse load, but only on the value of the integral Ι.


PERIODIC EXCITATION A SDF system is excited by a periodic (but not harmonic) load.

Newton’s equation is

u(t )

∞

k p(t )

m

po

j =1

mu&& + ku = ao → u p = To

2To

t

3To

∞

[

p(t ) = ao + ∑ a j cos( jωot ) + b j sin( jωot ) j =1

2π To

aj = To

1 ∫ p(t ) d t To 0

bj =

mu&& + ku = b j sin( jωot ) → u p =

]

To

2 ∫ p(t ) cos( jωot ) d t To 0 To

ao k

mu&& + ku = b j sin( jωot ) → u p =

Idea : a periodic function can be separated into its harmonic components using Fourier series.

ao =

]

The steady state response is calculated by using the theorem of superposition.

p(t )

ωo =

[

mu&& + ku = p(t ) = ao + ∑ a j cos( jωot ) + b j sin( jωot )

2 ∫ p(t ) sin( jωot ) d t To 0

b j /k 2

sin( jωot )

2

sin( jωot )

1 − ( jω o / ω n ) b j /k 1 − ( jω o / ω n )

The total steady state response is then

1/ k ao ∞ u(t ) = +∑ a j cos( jωot ) + b j sin( jωot ) k j =1 1 − ( jωo / ωn )2

[

]


Example

0 ≤ t ≤ To

p(t ) =

po t To

1 To po p ao = t dt = o ∫ To 0 To 2 2 To po aj = t cos( jωot ) d t = 0 ∫ To 0 To 2 To po p bj = t sin( jωοt ) d t = − o ∫ To 0 To jπ The total steady state response is then

po ∞ po / jπk u(t ) = −∑ sin( jωot ) 2k j =1 1 − ( jωo / ωn )2

with

ωo =2 ωn

u(t ) = +

po po  sin(ωot ) sin(2ωot ) − + 2k πk  − 3 − 30 sin(3ωot ) sin(4ωot )  + ... − 105 − 252 

Only the first three Fourier terms in u(t) must be considered to get an error less than 2 %. The plot of the load approximations with 2 and 3 Fourier terms shows that an inaccurate approximation of the load gives an accurate approximation of the response. The reason is that the higher frequencies in the load do not give any contributions to the response.


NON PERIODIC EXCITATION p(t ) u (t )

p(t )

t t

t To

Fourier series can also be used for non periodic loads. As example, the response u(t) of a damped SDF system loaded by the force p(t) is given above.

u(t )

t p(t) and therefore u(t) are non periodic. They become periodic by artificially adding the same model. Then a Fourier analysis can be done. The Fourier analysis of the artificial periodic problem is performed by only considering the steady state response. But the result obtained for the real problem includes both the transient and steady state response.

Remark: u(t) must be zero at the end of the period To.

Such an approach is not used for real problems, but it constitutes the basis for signal analysis, see the chapter experimental dynamics.


RAYLEIGH’S METHOD u E = max s E =0 k

Rayleigh’s method can be used to calculate approximately the lowest natural frequency of beams.

t

INTRODUCTION E =0 s E = max k

u(t ) Free vibrations of an undamped SDF system are considered.

u (t ) = uo sin(ω nt )

k m

Es max =

u& (t ) = uoω n sin(ω nt )

1 1 k u(t )2 = k uo2 sin2 (ωnt ) 2 2

kinetic energy

Ek (t ) =

Ek max =

1 m ωn2 uo2 2

conservation of energy

strain energy

Es (t ) =

1 k uo2 2

E = max s E =0 k

1 1 m u& (t )2 = m ωn2 uo2 cos2 (ωnt ) 2 2

Es max = Ek max

→

→

ωn =

1 1 k uo2 = m ωn2 uo2 2 2

k m


APPLICATION

strain energy

u( x, t )

M

EI L m

u ′′( x,t ) =

x

Es = Es beam = ∫

0

• using Es(t) = Ek(t) to get the natural frequency

Let’s take

 πx  u( x,t ) = Y sin  sin(ωt )  L 

which is the exact eigenmode if M = 0.

π2

 πx  Y sin   sin(ωt ) 2 L  L 

1 EI [u ′′( x,t )] 2 d x 2

4 L 1  πx  2 π 2 = EI Y 4 sin (ωt ) ∫ sin2   d x 2 L  L  0

The method consists in

• calculate Es(t) and Ek(t)

∂x 2

= −

L

The same approach can be used to calculate an approximative value of the lowest natural frequency of the beam above.

• estimating the vibration shape (eigenmode)

∂ 2u( x,t )

L

2  πx 

1L ∫ sin  L  d x = 2 ∫   0 0

L   2 πx   1 − cos d x =     2  L  

cos(2a) = 1 − 2 sin2 (a)

Es max

4 1 2 π = EI Y 2 2L3


kinetic energy

conservation of energy

∂u( x,t )  πx  u& ( x,t ) = = Y ω sin  cos(ωt ) ∂t  L 

Ek max = Es max

→

Ek = Ek M + Ek beam

Ek M =

1 1 2 M [u& (L / 2,t )] = M Y 2ω 2 cos2 (ωt ) 2 2

Ek beam

1  πx  m Y 2ω 2 cos2 (ωt ) ∫ sin2   d x 2  L  0

Ek max =

L

1 2 2 mL  Y ω M +  2 2  

EI L3(2M + mL )

particular cases

1 = ∫ m [u& ( x,t )] 2 d x 0 2 =

ω = π2

M=0 ω = π

L

L π4  → ω  M + m  = EI 3 2 2L  2

m=0

2

exact solution since the eigenmode is exact.

EI mL4

π2 ω = 2

EI ML3

= 6.98

EI ML3

In that case, the beam is considered as a spring with

k =

which gives as exact solution

48EI

ω =

k = M

ML3

48EI L3

= 6.93

EI ML3


REMARKS • Only an estimate of the natural frequency can be calculated. • The accuracy of the result depends entirely on the shape function which is assumed to represent the eigenmode.

a possibility is to take the deflected shape corresponding to the weight of the structure. (implemented in some f.e.m. codes) example M

F = Mg

• The natural frequency calculated by Rayleigh’s method is always greater than the exact value.

SELECTION OF THE SHAPE FUNCTION The shape function (eigenmode) must be kinematically admissible, i.e. must satisfy the displacements boundary conditions at the supports.

EI L m

q = mg

CONCLUSION The main interest of the Rayleigh’s method lies in its ability to provide useful estimation of the natural frequency from any reasonable assumption of the eigenmode.


SDF - TIME INTEGRATION METHODS u(t )

p

pi +1

p(t )

c k

pi

m

The response at time step i +1 is calculated from the equation of motion, a difference expression, and known responses at one or more preceding time steps. Equations of motion at time i and i +1.

t i t i +1

t

m u&&i + c u&i + k ui = pi

m u&& + c u& + k u = p(t )

The load p is time discretised.

p = [ p0 p1 p2 ... pi −1 pi pi +1 ... pn

]

The purpose is to calculate u (and u& , u&& if required) at the discrete time instants.

u = [ u0 u1 u2 ... ui −1 ui ui +1 ... un

]

u& = [ u&0 u&1 u& 2 ... u&i −1 u&i u&i +1 ... u&n

]

u&& = [ u&&0 u&&1 u&&2 ... u&&i −1 u&&i u&&i +1 ... u&&n

]

(1)

m u&&i +1 + c u&i +1 + k ui +1 = pi +1

(2)

An explicit algorithm uses a difference expression of the general form

ui +1 = f ( ui , u&i , u&&i , ui −1, u&i −1 ... ) which is combined with equation (1) An implicit algorithm uses a difference expression of the general form

ui +1 = f ( u&i +1, u&&i +1 , ui , u&i , u&&i , ui −1 ... ) which is combined with equation (2)


CENTRAL DIFFERENCE METHOD Introducing the two difference expressions in the equilibrium equation at time i

(explicit algorithm)

m u&&i + c u&i + k ui = pí

difference expressions :

u −u u&i = i +1 i −1 2 ∆t

u&&i =

ui +1 − 2 ui + ui −1

gives the unknown ui +1 as

( ∆t )2

kˆ ui +1 = pˆ i u

ui +1 ui ui −1

∆t

assumption u = at 2 + bt + c

∆t t

ti −1

ti

ti +1

m c kˆ = + ( ∆t )2 2 ∆t

 m  c  2m  pˆ i = pi −  − u k u − −  i −1  2 2 i 2 ∆ t ( ∆t )   ( ∆t )  

initialisation The given initial conditions are u0 and u&0 .

The displacement curve is assumed to be a parabola between 3 consecutive points.

But u0 and u−1 are required to determine u1.


The difference expressions for i = 0 are

u&0 =

u1 − u−1 2 ∆t

u&&0 =

Algorithm – central difference method

u1 − 2 u0 + u−1 ( ∆t ) 2

which after substitution gives

(∆t )2 u−1 = u0 − ∆t u&0 + u&&0 2

1 initial calculations

p − c u&0 − k u0 u&&0 = 0 m

(∆t )2 u−1 = u0 − ∆t u&0 + u&&0 2

m c kˆ = + ( ∆t ) 2 2 ∆t

m

a=

( ∆t )2

−

c 2 ∆t

b=k−

2m ( ∆ t )2

2 calculations for step i The equation of motion at time t = 0 gives the initial acceleration as

p − c u&0 − k u0 u&&0 = 0 m

pˆ i = pi − a ui −1 − b ui if required u&i =

ui +1 = pˆ i kˆ

ui +1 − ui −1 && ui +1 − 2 ui + ui −1 ui = 2 ∆t ( ∆t )2

3 repetition for the next time step replace i by i + 1 and repeat step 2


AVERAGE ACCELERATION METHOD (implicit algorithm)

Acceleration (constant)

difference expressions :

u&i +1 = u&i +

∆t (u&&i + u&&i +1 ) 2

u&&(τ ) =

1 (u&&i + u&&i +1 ) 2

u&&i +1

u&&i

(3)

( ∆t )2 ui +1 = ui + u&i ∆t + (u&&i + u&&i +1 ) 4

(4)

Velocity (linear)

u& (τ ) = u&i +

τ && && (ui + ui +1 ) 2

u&i +1 u&i

The acceleration is assumed constant between ti and ti+1. it is given by

1 u&&(τ ) = (u&&i + u&&i +1 ) 2

(3)

2 τ & && && Displacement u(τ ) = ui + ui τ + 4 (ui + ui +1 )

(Quadratic)

0 ≤ τ ≤ ∆t

ui +1 ui

The velocity is then linear and the displacement quadratic, which gives equations (3) and (4).

(4)

τ ∆t


The system constituted by equations (2) (3) (4) has 3 unknowns, u&&i +1 , u&i +1 and ui +1 . The solution of the system can be calculated as

kˆ ui +1 = pˆ i

2c 4m kˆ = k + + ∆t ( ∆t ) 2

 4u  4u&  2u  pˆ = pi +1 + c  i + u&i  + m  i 2 + i + u&&i  ∆t  ∆t   ( ∆t ) 

u&i +1 = u&&i +1 =

2 (ui +1 − ui ) − u&i ∆t

Algorithm – average acceleration method 1 initial calculations

u&&0 =

p0 − c u&0 − k u0 m

2 calculations for step i

 4u  4u&  2u  pˆ = pi +1 + c  i + u&i  + m  i 2 + i + u&&i  ∆t  ∆t   ( ∆t ) 

kˆ ui +1 = pˆ i

1 ( pi +1 − c u&i +1 − k ui +1 ) m

2c 4m kˆ = k + + ∆t ( ∆t ) 2

u&&i +1 =

u&i +1 =

2 (ui +1 − ui ) − u&i ∆t

1 ( pi +1 − c u&i +1 − k ui +1 ) m

The given initial conditions are u0 and u&0 . u&&0 can be calculated by

p − c u&0 − k u0 u&&0 = 0 m

3 repetition for the next time step replace i by i + 1 and repeat step 2


STABILITY 1.5

The central difference method will blow up, giving meaningless results, if the time step chosen is not short enough. The specific requirements for stability is :

1

0.5

0

-0.5

∆t 1 < Tn π

-1

The average acceleration method is unconditionally stable, which means that the procedure leads to bounded solutions regardless of the time step length.

-1.5 0

0.5

1

1.5

2

2.5

3

3.5

4

10 8 6 4

example

∆t = 0.01 ∆t = 0.07

∆t = 0.01 ∆t = 0.08

2

m u&& + c u& + k u = sin(2π t )

ωn = 8π

k =1

ξ = 0.05

0 -2 -4 -6

central difference method :

∆t < 0.0796

-8 -10

0

0.5

1

1.5

2


HOW TO CHOOSE ∆t ? ∆t must be small enough to get a good accuracy, and long enough to be computationally efficient.

For SDF systems, stability criteria are not restrictive since ∆t must be considerably smaller than the stability limit to ensure adequate accuracy. Stability of the numerical method is however important for the analysis of multiple degrees of freedom systems where it is often necessary to use unconditionally stable methods. One useful technique for selecting the time step is to solve the problem with a time step that seems reasonable, then repeat the solution with a smaller time step and compare the results, continuing the process until to successive solutions are close enough.

NONLINEAR RESPONSE If the structure undergoes large displacements or plastic deformation, the equation becomes

m u&& + c u& + f (u ) = p(t ) where the internal force f(u) is a nonlinear function.

The equation for the central difference method can be solved explicitly.

 m c  2m − kˆi ui +1 = pi − fi (ui ) −  u + u  i −1 2 2 i ∆ 2 t ∆ t t ( ) ( ∆ )   The equation for the average acceleration method must be solved by iteration (implicitly) since the first term in [] is unknown.

 ∂f 2c 4m  + + ui +1 = pˆ i +1  ∂u  2  i +1 ∆t ( ∆t ) 


2D CONTINUOUS BEAMS Free vibration in bending of a beam with distributed mass and flexibility is studied.

u( x,t ) ?

p( x,t ) x x

dx

EI m (kg/m)

M

V

∂t 2 dx

∂V V+ dx ∂x

∂M M+ dx ∂x

∂V   = V − V + x d   + pdx 2 ∂ x ∂t  

∂ 2u

m

∂ 2u ∂t 2

+

∂V =p ∂x

(1)

Rotational static equation ( the inertial moment associated with the angular acceleration is neglected ) : ∂M V= ∂x

M = EI

p dx

∂ 2u

m dx →

DIFFERENTIAL EQUATION

u

Vertical dynamic equation

∂ 2u ∂x 2

∂V ∂ 2M ∂ 4u = = EI 4 ∂x ∂x 2 ∂x

→

Equations (1) and (2) give

EI

∂ 4u ∂x

4

+ m

∂ 2u ∂t

2

= p(t )

(2)


FREE VIBRATION

m f&&(t ) φ (4) ( x ) → = − φ (x) EI f (t )

The differential equation

EI

∂ 4u ∂x

4

+ m

∂ 2u ∂t

2

The above equation must hold for every x and t

= 0

can be solved analytically.

φ (4) ( x ) m f&&(t ) → = − = a4 φ ( x) EI f (t )

→

A possible solution is of the form

u( x,t ) = φ ( x ) ⋅ f (t )

 φ (4) ( x ) − a 4φ ( x ) = 0    f&&(t ) + ω 2 f (t ) = 0 

(a constant)

EI a 4 ω = m 2

The solution is which gives

EI φ

with

φ

(4)

(4)

( x ) f (t ) + m φ ( x ) f&&(t ) = 0

(x) =

d4 φ ( x ) dx

4

f&&(t ) =

d2 f (t ) dt 2

 f (t ) = Ecos(ωt ) + Fsin(ωt )   φ ( x ) = Asin( ax ) + Bcos(ax ) + Csinh( ax ) + Dcosh(ax ) E, F are determined by the initial conditions. A, B, C, D are determined by the boundary conditions.


EXAMPLE: SIMPLY SUPPORTED BEAM  φ (0) = 0   φ ′′(0) = 0

u( x,t ) = φ ( x ) ⋅ f (t )

u

→

 B+ D = 0  − B+ D = 0

→

 B = 0   D = 0

EI L m x

Boundary conditions

M = EI

      

u( x = 0,t ) = 0 u( x = L,t ) = 0 M ( x = 0,t ) = 0 M ( x = L,t ) = 0

→

      

∂ 2u ∂x

2

 φ (L) = 0  A sin(aL) + C sinh( aL) = 0 →    φ ′′(L) = 0  − A sin(aL) + C sinh( aL) = 0  A sin(aL) = 0 →   C sinh( aL) = 0

= EI φ ′′( x ) ⋅ f (t )

φ (0) = 0 φ (L ) = 0

(2)

(2) → C = 0 or a = 0 (no motion) → C = 0

(1) → A = 0 (no motion) or sin(aL) = 0

φ ′′(0) = 0 φ ′′(L) = 0

(1)

Conclusion: the boundary conditions require

φ ( x) = Asin(ax) + Bcos(ax) + Csinh(ax) + Dcosh(ax) φ ′′( x ) = a 2 [− Asin( ax ) − Bcos(ax ) + Csinh( ax ) + Dcosh(ax )]

sin(aL) = 0 and B = C = D = 0

sin(aL) = 0 EIa 4 ω = m 2

B=C=D=0

→ →

→

aL = nπ

nπ a= L

→

ω n = n 2 π2


visualisation of the three first eigenmodes

EI mL4 nπx    L 

φn ( x ) = A sin

ωn are the natural circular frequencies φn are the natural eigenmodes An infinite number of solutions un have been found

un ( x,t ) = φn ( x ) ⋅ [ En cos(ωnt ) + Fn sin(ωnt ) ]

n=1

ω1 = π2

EI mL4

πx 

φ1( x ) = A sin

 L   n=2

EI

ω 2 = 4 π2

mL4  2πx    L 

φ2 ( x ) = A sin

n=3 The general solution for free vibration is a superposition of all the eigenmodes. ∞

u( x,t ) = ∑ φn ( x ) ⋅ [ En cos(ωnt ) + Fn sin(ωnt ) ] n =1

En Fn are given by the init . cond. u( x,t = 0) , u& ( x,t = 0)

ω 3 = 9 π2

EI mL4  3πx    L 

φ3( x ) = A sin


MDF – FREE VIBRATION Free vibration of multi degree of freedom systems are studied.

The coefficients in parenthesis are determined by considering first u1 alone and then u2 alone.

EXAMPLE 1 - TRAIN u1 k

→

 m u&&1 + 2k u1 − k u2 = 0   m u&&2 − k u1 + 2k u2 = 0

u2

k m

u1 > 0

k m

ku1

ku1

u2 = 0

ku1

Differential equations Newton’s law is applied for each coach

 m u&&1 = ( − 2k ) u1 + ( k ) u2   m u&&2 = ( k ) u1 + ( − 2k ) u2

u1 = 0

u2 > 0

ku2

ku2 ku2


Mass and Stiffness matrices The differential equations can be written as

0  m 0  u&&1   2k − k   u1  + =        0 m  u&&  − k 2k  u 0     2    2

mass matrix

stiffness matrix

[m ]

[k]

 2k  − k   

first column of [k]

u1 = 1

u2 = 0

2k

k

− k  second column of [k]    2k  The stiffness matrix can be determined by the classical way. As example, the components of first column of [k] are the external forces that must be applied to the masses in order to keep the system in equilibrium when a unit displacement u1=0 is applied.

u1 = 0

u2 = 1

k

2k


Eigenvalue problem 0  m 0  u&&1   2k − k   u1  + =    0 m  u&&   − k 2k  u  0     2    2

φ  φ  0 − ωn2 [ m]  1  sin(ωnt ) + [ k ]  1  sin(ωnt ) =   0 φ2  φ2  →

( [ k] − ωn2 [ m] ) φφ1  sin(ωnt ) = 00  

 2

The solution of this system is of the form This equation is valid for every t, which implies

u1  φ1  =     sin(ωnt ) u  2 φ2 

φ1    ? φ2 

ωn ?

Since φ1 , φ2 , ωn are constants, derivation gives

u&&1  2  &&  = − ωn u2 

φ1    sin(ωnt ) φ2 

By introducing the two expressions above in the system, it is obtained

( [ k] − ωn2 [ m] ) φφ1  = 00  2

 

(1)

eigenvalue problem

This system has one trivial solution (φ1 = φ2 = 0) which corresponds to equilibrium (no motion). Other solutions can be found if the following condition is respected:

det

( [ k] − ωn2 [ m] ) = 0

(2)


( 2)

→

2k − ωn2m −k = 0 2 −k 2k − ω n m

ω1 =

(2k − ωn2m)2 − k 2 = 0

→

(2k − ωn2m − k ) (2k − ωn2m + k ) = 0

ωn =

k m

or

ωn =

3k m

Two values of ωn have been found. For each of them, φ1 φ2 can be determined by introducing the value of ωn in the system (1).

( [ k] − ωn2 [ m] ) φφ1  = 00  2

(ω12m = k )

 k (1) →  − k

→

→

k m

− k k 

φ1  0    =   0  φ2 

This system has an infinite number of solutions, one of them is φ1   1 =      1 φ2 

ω2 =

3k m

(ω 22m = 3k )

− k (1) →  − k

− k − k 

φ1  0 =     0 φ2 

(1)

  → e . g.

φ1   1 =     − 1 φ2 


Conclusion : Two solutions have been found

{φ1}

ω1

u1(t )   1 =     sin( 3k m ⋅ t ) u ( t ) − 1  2 

ω1 ω 2

{ φ1} { φ2 }

 1   sin(ω2t + θ 2 ) − 1

Physical interpretation of the eigenmodes

u1(t )  1 =     sin( k m ⋅ t ) u t ( ) 1  2 

{φ2 }

 u1(t )  1 = C   1   sin(ω1t + θ1) + C2 u ( t ) 1  2 

Free vibration is initiated by an initial deflection corresponding to eigenmode 1:

u1(0) = A

u2 (0) = A

u&1(0) = 0

u& 2 (0) = 0

ω2

: natural circular frequencies : eigenmodes

Every linear combination of this two solutions is also a solution of the free vibration problem. The general solution can be written as

With these initial conditions, the response is

u1(t )  1 A =     sin(ω1t ) u ( t ) 1  2  • The motion of both solids is harmonic with ω1 for circular frequency. • The deflected shape is constant in time and corresponds to eigenmode 1.


EXAMPLE 2 – THREE STORY FRAME m

u1 = 1

L

L = 200 mm m

u2

first column

 k  − k   0

k

m = 2 kg

u1

k

    

L

columns

m

b = 10 mm

u3

k

L

h = 1 mm

E = 200 GPa

u2 = 1

second column 2k

k

( 4 columns by story )

k=

−k   2k −k 

    

bh3 I= 12

48EI L3

third column

m [ m] =  0  0

0 m 0

0 0  m 

 k [ k ] = − k  0

−k 2k −k

0 − k  2k 

( [k] is symmetric )

k

u3 = 1

2k

 0     −k   2k   


Eigenvalue problem

( [ k] − ω2 [ m] ) { φ } = {0} This is solved numerically by MATLAB.

M = [ m]

The eigenmodes are determined with an arbitrary multiplicative constant. By convention in dynamics, the eigenmodes as given on the form

φ11    φ1 = φ21    φ31 

with φ11 = 1

[F, E ] = eig(K ,M )

K = [ k]

The eigenmodes in MATLAB are calculated such

solves the eigenvalue problem

[ k] {φ } = λ [ m] {φ }

→ ω2 = λ The MATLAB solution is given as

ω12 0 0   E =  0 ω22 0  0 0 ω32   

 F =  

        φ1  φ2  φ3          

φ11    φ1 = φ21    φ31 

2 2 2 with φ11 + φ21 + φ31 = 1

In order to obtain the dynamic convention, the following scale command must be used.

F = [F (:,1) / F (1,1) F (:,2) / F (1,2) F (:,3) / F (1,3) ]


Results given by MATLAB

f1 = 1.58 Hz

f2 = 4.44 Hz

1.000   φ1 = 0.802    0.445

 1.000   φ2 = − 0.555    − 1.247

1.000

0.802

0.445

ω3 = 40.3 rad/s

ω 2 = 27.9 rad/s

ω1 = 9.95 rad/s

1.000

− 0.555

− 1.247

f3 = 6.41 Hz  1.000   φ3 = − 2.247    1.802 1.000

− 2.247

1.802

Physical interpretation of the eigenmodes If free vibration is initiated by initial displacements corresponding to the eigenmode n, the vibration of each story will be harmonic with frequency fn and the structure will vibrate with a constant deflected shape corresponding to the eigenmode n


Comparison between discrete and distributed systems

discrete systems eigenmode = vector number of natural frequencies = number of degrees of freedom

distributed systems

infinite number of natural frequencies

φ11    φ1 = φ21  φ   31 

eigenmode = function

πx φ1( x ) = A sin   L 

the eigenmodes are defined with a multiplicative constant

convention

φ11 = 1

convention

A = 1


MDF – MODE ANALYSIS Summary of the eigenvalue problem The mode superposition approach is used in order to calculate the response of a MDF system to an applied load. First harmonic loads and then arbitrary loads are considered. Finally damping is introduced.

[ m] =

m 0   0 m  

[ k] =

eigenvalue problem

 2k  − k

− k 2k 

( [ k] − ωn2 [ m] ) {φ} = {0}

HARMONIC LOAD – TRAIN The first coach is excited by a harmonic load

u1

u2

po sin(ωt ) m

k

m

k

k

ω1 = k m

{φ1} =   1 1

ω 2 = 3k m

The matrix [φ φ] is defined as

[φ] = [ {φ1}{φ2 } ] = 

1 1  1 − 1

{φ2 } = 

1  − 1


Differential equations

[ m]

by derivation, it is obtained

u&&1   u1   po  [ ] k + =  &&      sin(ωt ) u u 0  2  2

 m u&&1 + 2k u1 − k u2 = po sin(ωt )   m u&&2 − k u1 + 2k u2 = 0 The two equations are coupled and it is hence not possible to solve directly the system.

Modal coordinates

[φ]

 q&&1(t )   &&  q2 (t )

Introducing these two expressions in the differential equations and multiplying on the left side by [φ φ]T gives

q  q&&1  [φ] [ m][φ]  &&  + [φ]T [ k ][φ]  1  = [φ]T {p}sin(ωt ) q2  q2  T

which can be rewritten as

The idea is to change coordinates (unknowns) and calculate first q1(t) and q2(t) which are defined by

u1(t )    = u ( t )  2 

u&&1(t )   &&  = u2 (t )

q (t ) [φ]  1  q2 (t )

q  q&&1  [ M]  &&  + [ K]  1  = q2  q2 

[ M] = [φ]T [ m] [φ] [ K] = [φ]T [ k ] [φ]

{P} sin(ωt )

{P} = [φ]T {p}


After calculations it is obtained

[ M] =

 2m 0     0 2m 

[ K] =

 2k 0     0 6k 

{P} =

 po     po 

2m 0   q&&1  2k 0   q1   po   0 2m q&&  +  0 6k  q  =  p  sin(ωt )    2    2  o

 2m q&&1 + 2k q1 = po sin(ωt )   2m q&&2 + 6k q2 = po sin(ωt ) Since Both [M] and [K] are diagonal matrices, the two equations are uncoupled and can be solved separately. Each equation represents a SDF system. The natural circular frequencies for the first and second equations are ω1 and ω2 which have been calculated by the eigenvalue problem.

Each equation are solved by considering only the steady state response

po 2k  q ( t ) = sin(ωt ) 1 2  1 − (ω ω1 )    po 6k q ( t ) = sin(ωt )  2 2 ω ω ( ) 1 −  2

ω1 = k m ω 2 = 3k m

u1 and u2 can now be calculated with

u1(t )    = u ( t )  2 

[φ]

 q1(t )    q ( t )  2 

  po 2k po 6k   u ( t ) = + sin(ωt )  1 2  1 − (ω ω )2 1 − (ω ω 2 )    1    po 2k po 6k   u ( t ) = − sin(ωt )  2 2 2  1 − (ω ω 2 )    1 − (ω ω1 )


Conclusion   po 2k po 6k   + sin(ωt )  u1(t ) =  2 2 ( ) ( ) ω ω ω ω − − 1 1   1 2     po 2k po 6k   − sin(ωt )  u2 (t ) =  2 2 ( ) ( ) 1 − 1 − ω ω ω ω   1 2 

After some time (steady state response), the structure vibrates with the same frequency as the force.

• the amplitudes of the vibrations become infinite (undamped case)

If a MDF structure is excited by a harmonic force whose the frequency is one of the natural frequencies of the structure , then after some while (steady state response) • the structure vibrates with the same frequency as the applied force

Resonance the load frequency ω is close to the first natural frequency of the structure ω1

ω → ω1

• the structure vibrates with a deflected shape corresponding to eigenmode 1.

u1  1 p / 2k o sin(ωt )   =  2 1 1 − (ω ω1 ) u2 

• the deflected shape is the eigenmode associated to the natural frequency • the amplitudes of the vibrations are large This phenomenon is called resonance


HARMONIC LOAD – 3 STORY FRAME Applied harmonic force A harmonic force is applied to the second story.

m

1sin(2πt )

u1

m

u2

m

u3

m 0 0  [ m] =  0 m 0   0 0 m 

1.000 1.000 1.000 [φ] = 0.802 − 0.555 − 2.247  0.445 − 1.247 1.802

{p} sin(2πt )

Mode superposition

u1(t )    u2 (t ) = u (t )   3 

0  k −k [ k ] = − k 2k − k   0 − k 2k 

Eigenvalue analysis

ω1   9.951     ω2  = 27.88     ω3  40.29

0    1  sin( 2πt ) = 0   

 q1(t )  [φ] q2 (t ) q (t )   3 

u&&1   u1   &&  [ m] u2  + [ k ] u2  =     u&&3  u 3 

→

{p} sin(2πt )

 q&&1   q1   &&  [ M] q2  + [ K] q2  = q&&  q   3  3

{P} sin(2πt )


0 0  3.682 [ M] = [φ]T [ m][φ] =  0 5.726 0  0 18.59  0

Each SDF system is solved by only considering the steady state response, which directly gives

q1(t ) =

0  364.7 0 [ K ] = [φ]T [ k ][φ] =  0 4452 0  0 30184  0  0.802 {P} = [φ]T {p} =  − 0.555   − 2.247

The uncoupled system is

 3.682 q&&1 + 364.7 q1 = 0.802 sin(2 πt )   5.726 q&&2 + 4452 q2 = − 0.555 sin(2 πt )  18.59 q&& + 30184 q = − 2.247 sin(2 πt ) 3 3 

q2 (t ) = q3(t ) =

0.802/ 364.7

1 − (2 π 9.95)

2

sin( 2πt ) = 3.656 10− 3 sin(2 πt )

− 0.555/ 4452

1 − (2 π 27.9 )

2

sin(2 πt ) = − 1.313 10− 4 sin(2 πt )

− 2.247/ 30184 1 − (2 π 40.3)

2

sin(2 πt ) = − 7.629 10−5 sin(2 πt )

Transformation to real coordinates

 3.656 10 3  u1(t )   q1(t )       4 [ ] [ ] u ( t ) = q ( t ) = − 1 . 313 10 φ φ  2   2    sin(2 πt ) u (t )  q (t )  − 7.629 10 5   3   3    3.449 10−3  u1(t )     −3  u ( t ) = 3 . 177 10  2    sin( 2πt ) u (t )  1.653 10−3   3   


ARBITRARY LOADING – TRAIN

q1 = [ • • • • • ]

1

p (N)

q1

0.8

q2 = [ • • • • • ]

0.6

p(t) is applied to the first solid

0.4

10

q2

0.2

t (s) 0

1

0 -0.2 -0.4

[ m]

u&&1  u1   p(t ) [ ] + = k  &&       0  u2  u2 

-0.6 -0.8 -1

0

0.2

0.4

0.6

0.8

1

1.2

1.4

1.6

1.8

2

t (s)

Transformation to modal coordinates

q&&1  q  [ M]   + [ K ]  1  = {P} q&&2  q2 

→

 2m q&&1 + 2k q1 = p(t )   2m q&&2 + 6k q2 = p(t )

Each SDF equation is solved by using a time integration method (lesson 4). For this particular case, the exact solution can be found (lesson 3).

Transformation to real coordinates

 u1(t )    = u ( t )  2 

[φ]

 q1(t )    q2 (t )

must be performed for each time step

u1(t )  φ11 φ12  • • • • •  ∗ ∗ ∗ ∗ ∗ = =         ∗ ∗ ∗ ∗ ∗ u2 (t ) φ21 φ22  • • • • • 


SYSTEMS WITH DAMPING The equation system of MDF systems with viscous damping is

[ m] {u&&} + [ c] {u& } + [ k ] {u} = {p(t )} The transformation to modal coordinates can be performed as before, which gives

{u} = [φ] {q} [ M] {q&&} + [ C] {q& } + [ K ] {q} = {P(t )} with

[ C] = [φ]T [ c][φ]

Problems : • how to define [φ φ] ? • how to obtain by experiments [C] or [c] ? • useful transformation only if [C] is diagonal

Solution : 1) Natural frequencies and modes are calculated by neglected damping. This is correct if the damping is small. 2) Transformation to normal coordinates : a damping coefficient ξn is introduced for each mode. ξn is determined by experiment. This gives n uncoupled equations on the form.

Mn q&&n + Cn q&n + K n qn = Pn (t )

ξn =

Cn 2 K nM n

3) Each SDF equation is solved by time integration method

→

{q(t )}

4) Transformation to real coordinates

{u} = [φ] {q}


DAMPING MATRIX Method two : Rayleigh damping For certain problems, it is better to solve directly the coupled differential equations

[ m] {u&&} + [ c] {u& } + [ k ] {u} = {p(t )} instead of using the modal transformation (see next section). The damping matrix [c] is then needed. [c] is calculated from the modal damping ratios ξn which are determined by experiments.

The damping matrix is taken as

[ c] =

It can easily be shown that the modal coordinates transformation leads to a diagonal modal damping matrix and that

ξn =

Method one

a0 [ m] + a1 [ k ]

a0 1 a + 1 ωn 2 ωn 2

[c] is calculated from the modal matrix [C]

Cn ξn = 2 K nM n

[ C] = [φ]T [ c] [φ]

→

→

C1  [ C ] =  C2   C3  

[ c] = ( [φ]T ) [ C][φ] −1

This method supposes that all the damping ratios ξn are known, which is usually not the case.

Hence, [c] is calculated using only two modal damping coefficients. The two modes with specified ξn should be chosen to ensure reasonable values for the other damping ratios. In practice, the lowest modes and the third or fourth lowest ones are used to determine a0 and a1.


EFFECT OF THE DAMPING ON HARMONIC VIBRATIONS – TRAIN The first solid is excited by a force po sin(ωt ) Transformation to modal coordinates

 2m q&&1 + C1 q&1 + 2k q1 = po sin(ωt )   2m q&&2 + C2 q& 2 + 6k q2 = po sin(ωt ) The steady state response for each equation is on the form

po / 2k

[1 − (ω ω ) ]

2 2

1

tan(θ1 ) =

+ [2ξ1(ω ω1 )]

2ξ1(ω ω1 ) 1 − (ω ω1)2

 G1 sin(ωt − α1 )   u1(t )   q1(t )    = [φ]   =   u ( t ) q ( t )  G2 sin(ωt − α 2 )   2   2 

Conclusion : • Both solids vibrate with the same frequency as the applied force, but with a certain difference of phase. It means that u1(t) and u2(t) don’t reach their maximal values at the same time.

 q1(t ) = D1 sin(ωt − θ1)   q2 (t ) = D2 sin(ωt − θ 2 )

D1 =

Transformation to real coordinates : after some work, the response can be written as

2

• The resonance properties enounced without damping are still valid.


MDF - TIME INTEGRATION METHODS The system of differential equations can also be solved by time stepping methods. These methods are similar to the ones for SDF systems (lesson 4) by replacing scalars by matrices and vectors.

{u}0

(explicit algorithm) This method is based on two difference expressions and the equilibrium at time i.

{u& }i = {u}i +1 − {u}i −1 {u&&}i = {u}i +1 − 2 {u}i2+ {u}i −1

[ m] {u&&} + [ c] {u& } + [ k ] {u} = {p(t )} Initial conditions

CENTRAL DIFFERENCE METHOD

2 ∆t

{u& }0

( ∆t )

[ m] {u&&}i + [ c] {u& }i + [ k ] {u}i

= {p}i

combination of these 3 eqs. gives

 [m] [c]  {u} = {p} −  [m] − [c]  {u} − [k ] − 2[m]  {u} +   i +1  i −1  i  i 2 2 2 2 ∆ 2 ∆ t t ( ∆ t ) ( ∆ t ) ( ∆ t )      

using the following approximation

1 [c] ( {u}i −1 − {u}i +1 ) ≈ 1 [c] ( {u}i −1 − {u}i ) 2 ∆t ∆t

it is obtained

[m] {u} = {p} −  [m] − [c]  {u} − [k] − 2[m] + [c] {u}   i −1   i i +1 i 2 2 2

( ∆t )

 ( ∆t )

2∆t 



( ∆t )

Then, if [m] is diagonal, the solution of simultaneous equations is not required.

∆t 


AVERAGE ACCELERATION METHOD

(implicit algorithm) For this method, a set of uncoupled equations cannot be obtained by simplification since the equations are on the form

 2[c] 4[m]  [ ] {u}i +1 = {pˆ }i k + +  2 ∆ t ( ∆t )   And if [m] can often be taken as diagonal, this is not possible for the stiffness matrix [k]

CONCLUDING REMARKS The modal superposition requires that the system is linear, i.e. small elastic deformations. By comparison to time integration methods, the main interest of the modal analysis is that in most of structural dynamics problems, only a few modal contributions need to be solved.

Example

[ m] {u&&} + [ c] {u& } + [ k ] {u} = {p(t )}

10000 d.o.f.

→ 10000 modal equations

Mn q&&n + Cn q&n + K n qn = Pn (t ) Only the 10 first equations, corresponding to the 10 lowest natural frequencies need to be solved. q11 …. q10000 can be neglected and don’t need to be calculated.

 u1(t )   ....   u (t )  10  u (t )  ≈  11  ....    u10000 (t ) 

{u(t )} ≈ {φ1}q1(t )

 q1(t )   ....   q (t ) [φ]  10  0  ....    0 

+ {φ2 }q2 (t ) ... + {φ10 }q10 (t )


Implicit time integration methods are suited when many modes would be needed in the modal method. These methods are often unconditionally stable and the time step ∆t is only limited by accuracy. Explicit time integration methods are often used in combination with diagonal mass matrices which allows a small cost per time step. These methods are often conditionally stable and the time step must be very small so that stability is ensured for all the modes: even if the response in the higher modes is insignificant, it will blow up if the stability requirements are not satisfied relative to these modes. Contrary to implicit methods, explicit methods do not require an iteration process for nonlinear problems.

INTERPRETATION OF THE MODAL SUPERPOSITION  u1(t )    u2 (t ) = u (t )   3 

 q1(t )  φ11  φ21  φ31  [φ] q2 (t ) = φ12  q1(t ) + φ22  q2 (t ) + φ32  q3(t ) q (t )  φ       3   13  φ23  φ33  eigenmode 1

=

q1(t )

+

eigenmode 2

q2 (t )

+

q3(t )

eigenmode 3


EARTHQUAKE ANALYSIS

The ground acceleration can be registered using accelerographs, see the example below.

An introduction to earthquake analysis is given. Analysis by response spectra is presented.

SDF SYSTEMS A SDF system is subjected to a ground motion ug(t). The deformation response u(t) is to be calculated.

m (u&&g + u&&) + c u& + k u = 0 u&& + 2ξ ωn u& + ωn2 u = −u&&g (t )

u(t ) m

ug (t )

u(t )

fs (t )

fs (t ) = k u(t ) = m ωn2 u(t ) = m A(t )

A(t ) = ω n2u(t ) ≠ u&&(t )

c k/2

EQUIVALENT STATIC FORCE

k/2

Pseudo accelerati on fs(t) is the force which must be applied statically in order to create a displacement u(t).


REPONSE SPECTRA A response spectrum is a plot of maximum response (e.g. displacement, velocity, acceleration) of SDF systems to a given ground acceleration versus systems parameters (Tn , ξ). A response spectrum is calculated numerically using time integration methods for many values of parameters (Tn , ξ). Example : Deformation response spectrum for El Centro earthquake

D = max u(t ) t


Deformation, pseudo-velocity and pseudoacceleration response spectra can be defined and ploted on the same graphs

Peak Deformatio n

D = max u (t )

Peak Pseudo − velocity

V = ωn D

Peak Pseudo − acceleration

A = ωn2 D

ωn : natural circular frequency of the SDF system.

COMBINED D-V-A SPECTRUM


EXAMPLE

L=10 m

A water tank is subjected to the El Centro earthquake. Calculate the maximum bending moment during the earthquake.

ωn =

m = 10000 kg k = 98.7 kN/m

ξ = 2%

k 2π = 3.14 rad/s → Tn = =2s m ωn

 D = 7.47 ⋅ 25.4 = 190 mm Spectrum →   A = 0.191 ⋅ 9.81 = 1.87 ms −2

( obs : A = ω n2 D )

fs = k ⋅ D = 18.7 kN

Mmax = 187 kNm When the equivalent static force has been determined, the internal forces and stresses can be determined using statics.


RESPONSE SPECTRUM CHARCTERISTICS General characteristics can be derived from the analysis of response spectra.

Tn = 2 π m k Tn < 0.03 s : rigid system no deformation u(t) ≈ 0 → D ≈ 0

Tn > 15 s : flexible system no total displacement u(t) = ug(t) → D = ugo

The spectrum can be divided in 3 period ranges :

Tn < 0.5 s

: acceleration sensitive region

0.5 < Tn < 3 s : velocity sensitive region Tn > 3 s

: displaceme nt sensitive region


ELASTIC DESIGN SPECTRUM

Problem: how to ensure that a structure will resist future earthquakes.

The elastic design spectrum is obtained from ground motions data recorded during past earthquakes at the site or in regions with near-similar conditions

EXAMPLE


MDF SHEAR BUILDING

m1

u1

 u&1   u1  0  m1 0 0   u&&1 + u&&g         0 m2 0  u&&2 + u&&g  + [ c] u& 2  + [ k ] u2  = 0  0 0 m  u&& + u&&  u&3  u3  0 3  2  g

m2 u2

m3 u3

 u&1   u1   m1  m1 0 0   u&&1         0 m2 0  u&&2  + [ c] u& 2  + [ k ] u2  = − m2  u&&g (t )   u&3  u3  m3   0 0 m3  u&&3  Eigenvalue analysis

ω1 ω 2 ω 3

φ11    φ1 = φ21  φ31 

φ12    φ2 = φ22  φ32 

φ13    φ3 = φ23  φ33 

φ11 φ12 φ13  [φ] = φ21 φ22 φ23  φ φ φ   31 32 33 

Modal transformation

 m1  M1   q&&1  C1   q&1  K1   q1        T   q&&2  +  C2  q& 2  +  K 2  q2  = − [φ] m2  u&&g (t ) M2    M3  q&&3  C3  q&3  K3  q3  m3    

Mn q&&n + Cn q& n + K n qn = − Pn u&&g (t ) q&&n + 2ξ n ωn q&n + ω n2 qn = −Γ n u&&g (t )

n = 1,2 ,3

Γn =

Pn Mn

m1  {P } = [φ] m2  m3  T


RESPONSE HISTORY ANALYSIS Solve numerically

( u&&g (t ) known )

q&&n + 2ξ ωn q&n + ω n2 qn = −Γ n u&&g (t )

Summation of the modal contributions

n = 1,2 ,3

→

qn (t ) n = 1,2,3

φ13  u1(t )   q1(t )  φ11  φ12            u2 (t ) = [φ] q2 (t ) = φ21  q1(t ) + φ22  q2 (t ) + φ23  q3(t ) q3(t )  u3(t )  φ31  φ32  φ33 

Equivalent static forces

 f1(t )   u1(t )   q1(t )        f2 (t ) = [ k ] u2 (t ) = [ k ][φ] q2 (t ) f3(t )  u3(t )  q3(t ) 

Observation : The contribution of the higher modes can often be neglected.

φ13   f1(t )  φ11  φ12          f2 (t ) = [ k ] φ21  q1(t ) + [ k ] φ22  q2 (t ) + [ k ] φ23  q3(t ) f3(t )  φ31  φ32  φ33  The internal forces and stresses can now be determined using statics.

f1(t )

u1(t )

f2 (t )

u2 (t )

f3 (t )

u3 (t )

Vb (t )

Example : base shear force static Vb (t ) = f1(t ) + f2 (t ) + f3 (t ) Vb max = max Vb (t )


RESPONSE SPECTRUM ANALYSIS

( u&&g (t ) unknown ) Summation of the modal contributions The design response spectrum is used to estimate the behaviour of the shear building.

u1(t )   q1(t )      u2 (t ) = [φ] q2 (t ) u3(t )  q3(t ) 

The modal transformation gave (page 7.7)

q&&n + 2ξ n ωn q& n + ω n2 qn = −Γ n u&&g (t ) (1)

n = 1,2 ,3

which can be compared to the equation for a SDF (page 7.1)

u&& + 2ξ n ω n u& + ω n2 u = −u&&g (t ) (2) Dn , the maximum displacement for equation (2) is found on the design response spectrum. The maximum displacement for equation (1) is

qnmax = Γ n Dn

but

 u1max  Γ D     1 1 u2 max  ≠ [φ] Γ 2 D2  (3) u3max   Γ 3 D3 

because D1, D2, D3 do not occur at the same time.

(3) over estimates the maximal real displacements.


The summation of modal contributions is performed in a different way. The idea is to calculate first maximum displacements and equivalent static forces for each mode.

f1(m1)

u1(m1)

f1(m2)

u1(m2)

f2(m1)

u2(m1)

f2(m 2)

u2(m 2)

f3(m1)

u3(m1)

f3(m 2)

u3(m 2)

MODE 1

MODE 2

u1(m 1)  Γ1D1  φ11        u2(m 1)  = [φ]  0  = φ21  Γ1D1 u3(m 1)   0  φ31   

 u1(m 2)   0  φ12        u2(m 2)  = [φ] Γ 2D2  = φ22  Γ 2D2 u3(m 2)   0  φ32   

 f1(m1)  Γ1D1  φ11        f2(m1)  = [ k ][φ]  0  = [ k ] φ21  Γ1D1 f3(m1)   0  φ31   

 f1(m 2)   0  φ12        f2(m 2)  = [ k ][φ] Γ 2D2  = [ k ] φ22  Γ 2D2    0  φ32  f3(m 2) 

Same calculations for the third mode


Once a value of the required response has been calculated for each mode separately, the superposition of the different modes is performed using an approximative rule. The most common one is the SRSS approximation.

Maximum base shear force Maximum base shear force due to the mode 1

Vb(m1) = f1(m1) + f2(m1) + f3(m1) Maximum base shear force due to the mode 2

Maximum displacements

Vb(m 2) = f1(m 2) + f2(m 2) + f3(m 2)

 u2 + u12(m 2) + u12(m 3)  1 (m 1 )   u1max     2  2 2 u2 max  =  u2(m1) + u2(m 2) + u2(m 3)  u3max   2  2 2 u + u + u 3 (m 1 ) 3 (m 2 ) 3 (m 3 )  

Maximum base shear force due to the mode 3

Vb(m 2) = f1(m 3) + f2(m 3) + f3(m 3) Total maximum base shear force

Maximum shear force in the highest floor

Vb max = Vb2(m1) + Vb2(m 2) + Vb2(m 3)

f1

V1max = f12(m1) + f12(m 2) + f12(m 3)

V1

Observation : The contribution of the higher modes can often be neglected.


FINITE ELEMENT METHOD The equation system obtained by f.e.m. is

[ m] {u&&} + [ c] {u& } + [ k ] {u} = {p(t )}

The generalised load vector and stiffness matrix have been studied in the f.e.m. course and are identical in static and dynamic analyses. The damping matrix is not obtained by finite element discretisation. The different ways of introducing and dealing with damping have been studied in lesson 6.

with

[ m]

:

mass matrix

[ c]

:

damping matrix

[ k]

:

stiffness matrix

{p(t )}

:

generalise d load vector

{u(t )}

:

displaceme nt vector

The resolution of the equation system, by mode superposition or by direct time integration methods has been studied in lessons 4 and 6. The only remaining work is to present the derivation of the element mass matrix. The assembly of the element mass matrices to the structure mass matrix is performed in the same way as for the stiffness matrix.


CONSISTENT ELEMENT MASS MATRIX The equation of motion can be derived by equating the work done by the externally applied loads (external work) with the work absorbed by inertial and dissipative forces (internal work) for any virtual displacement (that is, for any imagined small motions that satisfies compatibility and essential boundary conditions.

Wi = ∫ v

{u} = [ N] {d} {u&&} = [ N] {d&&} {ε} = [ B] {d} Shape functions [N] are functions of space while nodal d.o.f. {d} are function of time.

T

{σ} = [ E] {ε} n

We = ∫ {δu} {F}dv + ∫ {δu} {Φ}ds + ∑ {δu}i { f }i T

v

Finite element discretisation provides

A linear elastic material is assumed

( {δε} {σ}+ {δu} ρ {u&&} ) dv T

{δu} and {δεε} represent virtual displacements and their corresponding strains.

T

s

T

i =1

where {F} and {Φ Φ} are the prescribed body forces and surface tractions, {f}i and {δu}i represent prescribed concentrated loads and their corresponding virtual displacements at a total of n points, ρ is the mass density.

The internal and external virtual works can then be rewritten as T  Wi = {δd}  ∫ v

We = {δd}

T

([ B] [ E][ B]{d}+ [ N] ρ [ N] {d&&} ) dv  T

T



n   T T T [ ] { } [ ] { } [ ] { } N F dv + N Φ ds + N f ∑ ∫ ∫ i i =1 s v 


The last equations in the preceding page are T  Wi = {δd}  ∫ v

We = {δd}

T

([ B] [ E][ B]{d}+ [ N] ρ [ N] {d&&} ) dv  T

This can be rewritten as

T

[ m] {u&&} + [ k ] {u} = {p(t )}



n   T T T [ ] { } [ ] { } [ ] { } N F dv N Φ ds N f + + ∑ ∫ ∫ i i =1 v s 

with

[ k ] = ∫ [ B]T [ E][ B] dv v

Wi = We is true for arbitrary {δd}, which gives

T &&}+ ∫ [ B]T [ E][ B] dv {d} = {p} ∫ ρ [ N] [ N] dv {d v

v n

{p} = ∫ [ N] {F}dv + ∫ [ N] {Φ}ds + ∑ [ N]T {f }i T

v

T

s

i =1

[ m] = ∫ ρ [ N]T [ N] dv v

The mass matrix defined in the last equation is called consistent element mass matrix. The word consistent emphasizes that this form has been derived using the same shape functions as the element stiffness matrix.


BERNOULLI 2D BEAM ELEMENT shape functions

v1

v2

θ1

θ2

EI m 1

L

v ( x,t ) =

[ N] {d}

2

mij = ∫ m ϕ i ( x ) ϕ j ( x ) dx 0

1

1 2

v1(t )  θ (t )  {d} =  1  v 2 (t ) θ 2 (t ) 

[ N] = [ ϕ1( x ) ϕ 2 ( x ) ϕ 3( x ) ϕ 4 ( x ) ]

L

m=ρA

x x ϕ1( x) = 1 − 3   + 2   L L

3

 x ϕ 2 ( x) = x 1 −   L

1 2

x x ϕ 3( x) = 3   − 2   L L

 156 22L 54 − 13L  2  13L − 3L2  mL  22L 4L  [ m] = 420  54 13L 156 − 22L   2 2   − 13L − 3L − 22L 4L 

3

2

1

x2  x  ϕ 4 ( x) =  − 1 L L 

 12 6L − 12 6L  2 2   EI  6L 4L − 6L 2L  [k] = 3 L  − 12 − 6L 12 − 6L   6L 2L2 − 6L 4L2   


BERNOULLI 2D BEAM ELEMENT WITH AXIAL DEFORMATIONS v1

θ1

v2

θ2

2

u2

EI A m 1

u1

L

A bar element is superposed to the previous beam element.

u( x,t ) = ϕ1( x ) u1 + ϕ 2 ( x ) u2 The shape functions for the bar element are

ϕ1( x ) = 1 −

0 70 0 0  140 0  0 156 22L 0 54 − 13L    2 2 0 22 L 4 L 0 13 L 3 L −   [ m] = mL   420 70 0 0 140 0 0    0 54 13L 0 156 − 22L   2 2   0 − 13L − 3L 0 − 22L 4L 

x L

ϕ2 (x) =

x L

− EA L 0 0 0 0  EA L   0 12EI L3 6EI L2 0 − 12EI L3 6EI L2    2 2 − 0 6 EI L 4 EI L 0 6 EI L 2 EI L   [ k] =   − EA L 0 0 EA L 0 0   3 2 3 2  0 − 12EI L − 6EI L 0 12EI L − 6EI L    6EI L2 2EI L 0 − 6EI L2 4EI L   0


LUMPED MASS MATRIX For certain applications it is better to use a diagonal (or lumped) mass matrix. Two obvious advantages are less storage space and less processing time, especially in case of an explicit time integration scheme.

The rotational inertia is defined by considering that a uniform slender bar of length L/2 and mass mL/2 is attached at each node. The associate inertia moment is

J = (mL / 2)(L / 2)2 / 3

PARTICLE MASS LUMPING Example : Bernoulli 2D beam element

v1

v2

θ1

θ2 m=ρA

EI m 1

L

2

An easy way of obtaining a lumped mass matrix is to replace the distributed mass by two particles of mass mL/2 at each node.

θ1

mL 2 1

u1

mL3 24

θ2

mL 2

mL3 24

0 1 2  mL 0 L / 12 [ m] = 2 0 0 0 0 

2

    0 L2 / 12

0 0 1

0 0 0

u2


HRZ LUMPING Different methods can be used to transform the consistent element mass matrix and obtain a diagonal matrix. One of them is the HRZ lumping. The idea is to compute only diagonal terms of the consistent element mass matrix, then scale them so as to preserve the total element mass. Since there may be both translational and rotational d.o.f., the method for an element of total mass M is 1. Compute diagonal coefficients mij of the consistent element mass matrix. 2. For each coordinate direction in which motion is described by the element d.o.f. a. determine a number S by adding the mii associated with translational d.o.f. only. b. multiply all coefficients mij associated with this direction by the ratio M/S.

Example : Bernoulli 2D beam element 1.

0 0  156 0 2  0 0  mL  0 4L  [ m] = 420  0 0 156 0   0 0 0 4L2  

2.a.

2.b.

312 S= mL 420 M = mL

0 1 2  mL 0 L / 39 [ m] = 2 0 0  0 0

    0 L2 / 39

0 0 1

0 0 0


EXAMPLE Mass and stiffness matrices for the structure The natural frequencies and eigenmodes of a simply supported beam are studied. The beam is discretised in 2 elements.

EIALm

m=ρA

θ1

θ2

3

2

v2 L/2

6Le − 12 6Le 0 0   12   2 2 6 L 4 L − 6 L 2 L 0 0 e e e  e    12+12 − 6Le+6Le − 12 6Le EI  − 12 − 6Le  [ k] = 3  2 2 2  2 Le 6Le 2Le − 6Le+6Le 4Le + 4Le − 6Le 2Le    0 0 − 12 − 6Le 12 − 6Le   2 2  0 6Le 2Le − 6Le 4Le   0

θ3

1

v1

Le = L / 2

v3 L/2

22 Le 54 − 13Le 0 0   156   2 2 4Le 13Le − 3Le 0 0   22Le   13Le 156+156 − 22 Le+22 Le 54 − 13Le mLe  54  [ m] = 2 2 2 420  − 13Le − 3L2e − 22 Le+22 Le 4Le + 4Le 13Le − 3Le     0 0 54 13Le 156 − 22 Le   2 2  0 0 − 13 L − 3 L − 22 L 4 L   e e e e 


Boundary conditions

• •  • •  •  •

v1 = v 3 = 0 v&&1 = v&&3 = 0 p1(t ) , p3(t ) unknowns

• • • • • • • • • •  • • • • • • • • • • • • • • •  • • • • •

v&&1  • θ&&  • 1    & & v 2  •  &&  +  θ 2  • v&&3  •  &&   θ •  3

θ1

θ2

1

Equation system

 4L2e  mL  13Le  420 − 3L2 e   0 

13Le

− 3L2e

312

0

0

8L2e

− 13Le

− 3L2e

  − 13Le   2  − 3Le  2  4Le  0

v1   p1(t )  θ   0  1     v 2   0    =   θ 0  2    p3(t ) v 3      θ 0    3

θ3 3

2

v1

• • • • • • • • • •  • • • • • • • • • • • • • • •  • • • • •

v2

 4L2e  θ&&1     v&&2  EI − 6Le  + 3  θ&&2  L  2L2e    & &  0 θ3  

v3

− 6Le

2L2e

24

0

0

8L2e

6Le

2L2e

0   6Le   2 2Le  2 4Le 

 θ1    v 2   = θ2    θ3 

0    0    0    0 


Solution of the eigenvalue problem

ω1 = 9.91

EI mL4

ω2 = 43.8

(9.87)  1  0.32 L    φ1 =   0    − 1 

(exact values in parenthesis)

EI mL4

(39.5)

mL4

θ 3 = −1

θ1 = 1

 1  − 0.05 L    φ3 =   0    − 1 

θ 2 = −1

ω4 = 201

EI mL4

(158)

eigenmode 2

θ2 = 0

v 2 = 0.32L

EI

(88.9)

 1  0   φ2 =   − 1   1 

eigenmode 1

θ1 = 1

ω3 = 110

 1 0   φ4 =    1  1 

eigenmode 3

θ3 = 1

θ1 = 1

θ2 = 0

v2 = 0

v 2 = −0.05L

θ 3 = −1


EXPERIMENTAL DYNAMICS The purpose of experimental dynamics is to determine by measurements the dynamics characteristics of a structure. This is performed by measuring the response of the structure to a known excitation. known excitation

known response unknown dynamics characteristics

Different types of excitations can be used

Free vibration : the structure is disturbed from equilibrium and then vibrates without any applied forces. Harmonic force : a sinusoidal force is applied through a special device. Pulse force

: a hammer is used to apply a short pulse force.

Unknown dynamics characteristics SDF : natural frequency

fn

Random force : random noise as e.g. traffic or wind loads are applied.

damping coefficient ξ MDF : natural frequencies

f1 f2 f3 …

eigenmodes

φ1 φ2 φ3 … ξ1 ξ2 ξ3 …

modes damping

For practical reasons, accelerations are easier to measure than displacements


SDF SYSTEMS

u(t )

Objective : determine experimentally fn and ξ

c k

m

u&&

C e −ξωnt u&&n u&&n + p

Free vibration test

t n+p

tn

u(t ) = C e−ξωnt sin(ω Dt + θ ) → u&&(t ) = D e

−ξωnt

sin(ωDt + α )

After derivation a similar expression for the acceleration is obtained and consequently the results derived in Lesson 1 for the displacements can be used.

TD =

2π

ωD

ξ < 0.1 →

ξ=

1 u&& ln n 2πp u&&n + p

→

fn =

1 1 ≈ Tn TD

t


u

Harmonic force test u max

The structure is excited by a harmonic force whose frequency is slowly increased step by step in order to reach the steady state response at each increment. The amplitudes of accelerations are measured at each step.

u max 2

f =

Results from Lesson 2 :

p(t ) = po sin(ωt ) u =

u = u sin(ωt + θ )

n

2 2

2 + [2ξ (ω ωn )]

ξ < 0.1 → fn = fr

ξ=

f2 − f1 f2 + f1

In reality, the maximal value is not known and a curve fitting is done.

By derivation, it is obtained :

u&& = − u&& sin(ωt + θ )

2π

f1 fr f2

po k

[1 − (ω ω ) ]

ω

u&& = ω 2 u

If the damping is law, the same method can be used with u&& instead of u .


Ambient vibration test The structure is excited by a small random force, which means that the force is not known. This method is often used for bridges where wind and traffic loads are used as random forces.

u&&(t ) = C e −ξωnt sin(ω Dt + α ) + Particular solution initial conditions u&&0 u&0

random load

The response (acceleration) is the sum of the free vibration solution and the particular solution which depends on the random load. The idea of this approach is to remove the particular solution and the initial condition u&0 by taking the average of many subrecords of same length ∆t and same u&&0 . The result is a free vibration solution which can be then used. The trigger value u&&0 and the length ∆t must be carefully chosen .

∆t

u&&

u&&0 t


MDF SYSTEMS Objective : determine experimentally f1 f2 f3 φ1 φ2 φ3 and ξ1 ξ2 ξ3

u1 u2 u3

Harmonic force test

ω = ω1 ( f = f1 )  u&&1(t )  φ11  A1 large     → → u&&2 (t ) = φ21  A1 sin(ωt + θ ) A 2 , A 3 negligible     u&&3(t )  φ31  The structure vibrates with a deflected shape corresponding to eigenmode 1 and the amplitudes of vibrations are large : resonance. u&&i

The structure is excited by a harmonic force whose frequency is slowly increased step by step in order to reach the steady state response at each increment. The amplitudes of accelerations are measured at each step. Results from Lesson 6 :

 u&&1(t )  φ11 φ12 φ13   A1        & & u2 (t ) = φ21 φ22 φ23   A 2  sin(ωt + θ )      & & φ31 φ32 φ33   A 3  u3(t ) 

u&&1 u&&2 u&&3

φ11    φ1 = φ21  φ   31 

φ11 φ21

φ31 f1

f2

f3

f =

ω 2π


The sign of φij is determined by the phase which is also measured.

Explication : the solution at the first resonance is e.g.

u&&1   1.00     o u&&2  = C − 0.55 sin(ωt + 30 )  &&    . u − 1 25  3  

 1.00 sin(ωt + 30o ) u&&1   which can be   & & u = C 0.55 sin(ωt + 30o + 180o )    2 rewritten as  &&   o o  1.25 sin(ωt + 30 + 180 ) u3 

Remark : this method assumes that only one mode is acting at each resonance. For that, the natural frequencies must be well separated:

ω = ωi → Ai large and other A j negligible Ai

A2

A1

A3

The measured phase is either 30° or 210° , which gives the sign of φij . f1

ξ1, ξ2, ξ3 are determined by band width method for each resonance.

f2

f

f3

In that case, the influence of the 2nd eigenmode can not be neglected at the 3rd resonance.


Free vibration test (damping neglected)

Results from Lesson 5 :

 u&&1(t )  φ11  φ12  φ13   &&        u2 (t ) = φ21  sin(ω1t + θ1 ) + φ22  sin(ω 2t + θ 2 ) + φ23  sin(ω3t + θ 3 ) u&& (t )    φ     3  φ31   32  φ33 

u&&1(t ) = φ11 sin(ω1t + θ1 ) + φ12 sin(ω 2t + θ 2 ) + φ13 sin(ω 3t + θ 3 ) u&&1

It is impossible to identify the dynamics parameters φ11 φ12 φ13 ω1 ω2 ω3 θ1 θ2 θ3 directly from the record of the acceleration of the first story. 0

The solution is to transform the record from the time domain to the frequency domain with a Fourier transformation:

t

∞

∞

n =1

n =1

u&&1(t ) = ∑ [an cos(ωnt ) + bn sin(ωnt )] = ∑ cn sin(ωnt + αn )


A numerical Fast Fourier transformation (FFT) gives the following two graphs, which are the representation of the signal in the frequency domain.

αn

cn 40

fn =

∞

u&&1(t ) = ∑ cn sin(ωnt + αn ) n =1

30

ωn

2π

30 o

10o

20

1.58

1.58

4.43

6.41

fn

6.41

- 30o

The required dynamics parameters can be easily determined form these two graphs:

u&&1(t) = 30 sin( 2π ⋅ 1.58t + 30o ) + 40 sin( 2 π ⋅ 4.43t − 30o ) + 20 sin( 2 π ⋅ 6.41t + 10o )

f1 = 1.58 Hz f2 = 4.43 Hz f3 = 6.41 Hz

φ11 = 30 φ12 = 40 φ13 = 20

θ1 = 30o θ 2 = − 30o θ 3 = 10o

fn


The eigenmodes can be entirely determined by doing the same operation for the 2nd and 3rd story.

 u&&1(t )  φ11  φ12  φ13   &&        u2 (t ) = φ21  sin(ω1t + θ1 ) + φ22  sin(ω 2t + θ 2 ) + φ23  sin(ω3t + θ 3 )   φ  φ  φ  u&&3(t )   31   32   33 

u&&1

u&&2

φ11 φ12 φ13

u&&3

φ32 φ33 φ31

φ23 φ21 φ22 f1

f2

f3

f

f1

f2

f3

f

f1

f2

f3

The sign of each φij is determined by the phase which is also given by a FFT.

f


Free vibration test (with damping) The same method is applied. The signal becomes periodic by artificially adding the same model (see Lesson 3) u&&1

u&&1

FFT

t

f f1

f2

f3

ξ1, ξ2, ξ3 are determined by band width method or by cutting the signal in the frequency domain and performing an inverse FFT to obtain free vibration of a SDF system corresponding to the eigenmode.


Pulse force excitation

p(t )

A hammer is used to apply a short pulse force. u&&1(t )

p(t )

t t

To

t

p(t) and the response are non periodic. They become periodic by artificially adding the same model. The artificial response is then the steady state response of the structure loaded by the artificial periodic load. The artificial load can be considered as the sum of harmonic terms. Each of these harmonic terms gives an harmonic response. ∞

part (t ) = ∑ pn sin(ωnt + α n )

u&&1(t )

t

ωn = ω1

D1n  φ11      → resonance D2n  = φ21      D3n  φ31 

n =1

u&&1(t )  ∞  &&  u2 (t ) = ∑ n =1 u&& (t )   3 art

D1n    D2n  sin(ωnt + θ1 ) D   3n 

The dynamic parameters can then be obtained by performing a FFT.


pn ∞

part (t ) = ∑ pn sin(ωnt + αn ) n =1

fn =

ωn

2π

The load obtained by the hammer must contain all the frequencies that have to be studied.

u&&i ( Din ) u&&1

u&&1(t )  D1n  ∞  &&    u2 (t ) = ∑ D2n  sin(ωnt + θ1) n =1   u&& (t )   3 art D3n 

u&&2 u&&3

φ11    φ1 = φ21  φ   31 

φ11 φ21

The FFT of the response is similar to the graph obtained with an harmonic force test (see page 9.5 )

φ31 f1

f2

f3

fn =

ωn

2π


SIGNAL ANALYSIS - FFT A signal can be represented in two different ways: in the time domain or in the frequency domain. Each domain has its own interest. The transformation between the two domains is performed numerically by a FFT or an inverse FFT. ∞

FFT

Y(f)

y(t)

 an (fn )   bn (fn )

or

y (t ) = ∑ [an cos(2 πfnt ) + bn sin(2 πfnt )]

 cn (fn )  αn (fn )

n =1 ∞

= ∑ cn sin( 2 πfnt + α n )

inv. FFT

n =1

The signal in the frequency domain is represented by two graphs (an bn) or (cn αn).

y

FFT

N points

t 0

T

αn

cn

N/2 points

N/2 points

inv. FFT 0

∆f

fmax

fn

0

∆f

fmax

fn


y

FFT

N points

t 0

0

y(t) is defined numerically by N points during a time T.

rate =

N T

m = 1,2,3…

N/2 points

N/2 points

inv. FFT

T

The FFT requires N = 2m

αn

cn

fmax

∆f

fn

0

fmax

∆f

cn and αn are defined by N / 2 points.

∆f =

1 T

fmax = ∆f

N rate = 2 2

In practice, ( ∆f , fmax ) are chosen and ( T , rate ) are then determined such that N = 2m.

fn

STRUCTURAL DYNAMICS – 5C1840

load F(t) • displacements ? structure

restriction : linear structures

• strains

?

• stresses

?

m u&& + c u& + k u = F (t )

2D continuum beams • Rayleigh’s method • pulse and periodic loads

SDF systems

F.E.M. • mass matrix

MDF systems

Earthquakes

• eigenvalue problem • mode superposition • free vibration • harmonic loads • time step methods

Experimental Methods

free vibration

eigenvalue problem

damping

eigenmodes

natural frequency

mode superposition

resonance

mode damping

transient response

Rayleigh damping

steady state response

consistent mass matrix

explicit time integration

lumped mass matrix

implicit time integration

time domain frequency domain

Structural Dynamics

Recommend Documents