Advanced Structural Dynamics and Vibration Control

Random Vibrations

Kishor Timsina, PhD

INTRODUCTION:

Let’s consider two major types of environmental loading on civil engineering structures:
wind loading and seismic loading.

Seismic loading is proportional to ground acceleration, which is random in nature:

Figure 1: North-south component of ground acceleration recorded at El Centro during the El Centro,
California earthquake of May 18, 1940 [Ref.16]
Wind Loading (drag force) is proportional to (wind velocity) 2, which is also random:

Figure 2: Wind speed at three heights on a tall mast in open terrain, East scale, Australia
(After Deacon, as quoted in ref.17)
If the time history of a loading is known, the structural response to the loading can be estimated from
the theory of structural dynamics.

However, it is normally impossible to know such future loading.

Then, how can we estimate the structural response?

How can we design our structures to withstand future unknown environmental loading?

One possible engineering solution:

From the available (past) records of loading, we extract a law(s) and/or parameters(s) which
characterize the loading.

These law(s) & parameter(s) are usually expressed in terms of statistical functions or parameters,
such as probability distributions, mean value, and variance.

In many cases , it is more physically meaningful to convert the loading in time domain into the loading
in frequency domain (see Figure 3)
Fourier Transformation Fourier Amplitude
Periodic function
Figure 3: Hypothetical examples of
zero-mean time record and
corresponding frequency spectra.

(a) Superposition of two sinusoids with

different frequencies f1 and f2 and
Random function Power Spectral Density
unequal amplitude. (narrow-band type)

(b) Superposition of several sinusoids

with respective frequencies clustered
within a narrow band around a
dominant frequency.

(c) Superposition of sinusoids with Random function Power Spectral Density

(wide-band type)
respective frequencies spread on a
rather wide band that includes a
“dominant” frequency.

Note that power spectral density can be seen as a probability distribution

function presented in the frequency domain.
The theory of random vibrations has been formulated basically in the frequency domain. The theory
provides a convenient framework for the analysis of the structural response to random loading.

RANDOM FORCE STRUCTURES RANDOM RESPONSE

Statistical model(s) mass, stiffness, damping Statistical model(s)
Statistical parameter(s) Statistical parameter(s)

These statistical models and parameters of response are very useful for the design of structures
since many physically meaningful response quantities can be derived.

These response quantities include: mean response, root mean square of response, peak response,
mean fatigue life, etc.

Physical insight can be gained. Computational requirement is minimum.

Probability Theory

7
Probability Theory

𝑃 is a probability function.
𝑃 is associated with an event.
𝑃(event) = 0 No chance that the event will become true.

𝑃(event) = 1 100% sure that the event will become true.

In general, 0≤P≤1
𝑃 (a major earthquake strikes Kathmandu in this year) = 0.001
𝑃 (the strength of a sample concrete cylinder exceeds 300 kg/cm2) = 0.30
𝑃 (the strength of a sample concrete cylinder exceeds 3000 kg/cm2) = 0 ?
𝑃 (the strength of a sample concrete cylinder exceeds 30 kg/cm2) = 0.999
𝑃 (you will get GPA 4.0 in this semester) = 0.4 ???
Let 𝑋 be a random variable. Let 𝑥 be a specified value (not random, it is a number that we specify
or that we are interested).

Cumulative Distribution Function (CDF) 𝐹 𝑥 ≡ 𝑃(𝑋 ≤ 𝑥)

The probability that a sample value of 𝑋 falls below the specified value 𝑥.

Probability Density Function (PDF) 𝐹 𝑥

𝑓 𝑥 𝑑𝑥 ≡ 𝑃 𝑎 < 𝑋 ≤ 𝑏

By these definitions, 𝐹 𝑥 = 𝑓 𝜉 𝑑𝜉

𝑓 (𝑥)
𝐹 ∞ = 𝑓 𝜉 𝑑𝜉 = 1

and, 𝑑𝐹 𝑥
= 𝑓 (𝑥)
𝑑𝑥
Suppose that we have 𝑁 sample values of 𝑋: 𝑥 , 𝑥 , 𝑥 , … … , 𝑥

1 →
Sample mean 𝑋= 𝑥 𝜇 The mean value of X
𝑁

1 →
Sample variance 𝑆 = 𝑥 −𝑋 𝜎 The variance of X
𝑁−1
Computed numbers from the sample data Actual statistical parameters of X

Ensemble average (or “expectation”) 𝐸[ ] a statistical operator

1 1
𝐸𝑋 = lim 𝑥 = 𝜇 = lim 𝑥 . = 𝑥 𝑓 𝑥 𝑑𝑥
→ 𝑁 → 𝑁
 1 1
𝐸𝑋 = lim 𝑥 = 𝜇 = lim 𝑥 . = 𝑥 𝑓 𝑥 𝑑𝑥
→ 𝑁 → 𝑁

Let g be a function of 𝑋 Examples: 𝑔(𝑋) = 𝑋 or 𝑔 𝑋 = sin(𝑋 )

1 1
𝐸 𝑔(𝑋) = lim 𝑔(𝑥 ) = lim 𝑔(𝑥 ) . = 𝑔 𝑥 𝑓 𝑥 𝑑𝑥
→ 𝑁 → 𝑁

𝐸𝑋 = 𝑥 𝑓 𝑥 𝑑𝑥 = The mean square of 𝑋

𝐸[𝑋 ] = the root−mean−square of 𝑋 (the rms value)

𝜎 ≡ 𝐸 𝑋−𝜇 = 𝑥−𝜇 𝑓 𝑥 𝑑𝑥 = 𝐸 𝑋 −𝜇
Guassian Random Distribution Many physical random quantities are known to be Gaussian type.

Let 𝑋 be a Gaussian random variable, that is

1 − 𝑥−𝜇
𝑓 (𝑥) is a Gaussian distribution. 𝑓 𝑥 = exp for − ∞ < 𝑥 < ∞
2𝜋. 𝜎 2𝜎

𝑓 (𝑥)
𝜎 𝜎 Symmetrical with respect to 𝜇 , bell shape, area under
the curve is “1”
𝜎 is a measure of scatterness around the mean 𝜇
(𝜎 is some time called “standard deviation”).

For a Gaussian random variable, the mean value 𝜇 and the standard deviation 𝜎 completely
describe its probability density function (2-parameter model).
 Random Process:

𝑋(𝑡) : a random function of time generated by a random process

Example: wind velocity fluctuation, temperature variation, etc.

The whole family or ensemble

of sample time histories of 𝑋 (𝑡)
which have been recorded.

.
.
.

A single individual time history 𝑥 ( ) (𝑡) called a sample function (see Fig. 4)

All sample functions (𝑥 (𝑡), where 𝑖 = 1, 2, 3, … ) are generated by the same physical process.
Figure 4: Schematic representation of a random
process 𝑋 𝑡 . Each 𝑥 ( ) (𝑡) is a sample function
of the ensemble.
Look at Figure 4. Let 𝑋 (𝑡 ) be a random variable 𝑋 .

A set of sample data for 𝑋(𝑡 ) is 𝑥 𝑡 , 𝑥 (𝑡 ), 𝑥 (𝑡 ) … … … … . . , 𝑥 (𝑡 )

This data set is obtained from 𝑁 records (𝑁 sample time histories).

1
Sample mean 𝑋(𝑡 ) = 𝑥 (𝑡 )
𝑁
1
(Actual) mean 𝜇 = 𝜇𝑡 = 𝐸 𝑋(𝑡 ) = lim 𝑥 (𝑡 )
→ 𝑁

Variance 𝜎 ( ) = 𝜎 𝑡 = 𝐸 𝑋(𝑡 ) − 𝜇(𝑡 )

1
Variance 𝜎 𝑡 = lim 𝑥 𝑡 − 𝑋(𝑡 )
→ 𝑁−1
Stationary Random Process:

Stationary means that the probability density function (and all other statistical functions) of the
process is time-invariant (i.e. not dependent on time).

1 𝑥−𝜇 𝑡
Non-stationary Gaussian process: 𝑓 𝑥 = exp −
2𝜋 𝜎(𝑡 2𝜎 𝑡

1 𝑥−𝜇
Stationary Gaussian process: 𝑓 𝑥 = exp −
2𝜋 𝜎 2𝜎
A universal distribution, independent of time

𝜇 𝑡 →𝜇
𝜎 𝑡 →𝜎

𝑓 𝑥 →𝑓 𝑥
Temporal average: All averages discussed so far have been ‘ensemble averages’.

Given a single sample record 𝑥 (𝑡) it is possible to obtain another kind of average called
‘temporal average’.
𝑥 (𝑡)

𝑡
𝑂
𝑇
() 1
The temporal mean 𝜇 ≡ 𝑥 𝑡 𝑑𝑡
𝑇

() 1
The temporal mean square 𝑋 ≡ {𝑥 ( ) (𝑡) 𝑑𝑡
𝑇

() 1
The temporal variance 𝜎 ≡ {𝑥 ( ) (𝑡) − 𝜇( ) 𝑑𝑡
𝑇
Ergodic Process

Within the subclass of stationary random processes there exists a further subclass
known as “ergodic process”.

An ergodic process is one for which ensemble averages are equal to the corresponding
temporal averages.

()
𝜇 = 𝜇
()
𝐸𝑋 = 𝑋 for 𝑇 → ∞

()
𝜎 = 𝜎

All important statistical information can be obtained from a single sufficiently long record.
It is common practice to assume that the process is ergodic.
Correlation between two random variables x and y
 1
Autocorrelation Function 𝑅 𝑡 ,𝑡 ≡ 𝐸 𝑋 𝑡 .𝑋 𝑡 = lim 𝑥 𝑡 .𝑥 𝑡
→ 𝑁
For special case where 𝑡 = 𝑡
For a zero-mean process,
that is, 𝜇 𝑡 = 0
𝑅 𝑡 ,𝑡 ≡ 𝐸𝑋 𝑡 = 𝜎 𝑡 +𝜇 𝑡 = 𝜎 (𝑡 )

If 𝑋(𝑡) is a stationary random process with zero mean, then

𝑅 𝑡 ,𝑡 = 𝐸 𝑋 𝑡 . 𝑋(𝑡 ) = 𝐸 𝑋 𝑡 . 𝑋(𝑡 − 𝑡 + 𝑡 ) let 𝜏 ≡ 𝑡 −𝑡

S
= 𝐸 𝑋 𝑡 . 𝑋 𝑡 +𝜏 = 𝐸 𝑋 𝑡 . 𝑋 𝑡 +𝜏 (the process is “stationary”)

𝑅 𝑡 ,𝑡 = 𝑅 𝑡 −𝑡 = 𝑅𝜏 A function of time lag 𝜏

𝑅 −𝜏 = 𝐸 𝑋 𝑡 . 𝑋 𝑡 − 𝜏 𝐸 𝑋 𝑡 +𝜏 . 𝑋 𝑡 = 𝐸 𝑋 𝑡 .𝑋 𝑡 + 𝜏 =𝑅 𝜏
𝑡 ≡𝑡−𝜏
𝑅0 = 𝜎 ≥ 𝑅𝜏 This can be proved by the probability theory
 𝑅𝜏
𝑅 𝑡 ,𝑡 ≡ 𝐸 𝑋 𝑡 .𝑋 𝑡 = 𝑅 𝑡 −𝑡 = 𝑅𝜏

𝑅 −𝜏 = 𝑅𝜏

𝑅0 = 𝜎 ≥ 𝑅𝜏

𝑋𝑡
For a stationary ergodic process,

/
1
𝑅 𝜏 = lim (𝑥 𝑡 . 𝑥 𝑡 + 𝜏 ) 𝑑𝑡 a time average
→ 𝑇
/

𝑥𝑡

𝑇/2 𝑇/2
Calculation of autocorrelation
Schematic diagram for measurement of autocorrelation function 𝑅(𝜏) for sample function 𝑓(𝑡).
The adjustable distance 𝑎 between tape heads is proportional to the lag 𝜏.
Schematic diagram of measurement of spectral density of sample function 𝑓(𝑡).
Fourier Analysis: Basic concept: a periodic function can be broken down to its harmonic
components.

𝑓 (𝑡) is a periodic function of time 𝑡 with period 𝑇,

that is, 𝑓(𝑡 + 𝑇) = 𝑓(𝑡) for any time 𝑡.

/
1
Where 𝑎 = 𝑓 𝑡 𝑑𝑡 (2a-c)
𝑇
/
/
The Fourier Series of 𝑓 𝑡 is given by 2 2𝑛𝜋
𝑎 = 𝑓 𝑡 . cos 𝑡 𝑑𝑡
𝑇 𝑇
2𝑛𝜋 2𝑛𝜋 /
𝑓(𝑡) = 𝑎 + 𝑎 cos 𝑡 + 𝑏 sin 𝑡
𝑇 𝑇 /
2 2𝑛𝜋
(1) 𝑏 = 𝑓 𝑡 . sin 𝑡 𝑑𝑡
𝑇 𝑇
/
Complex form of Fourier series:
The complex form is more compact for mathematical operations.

Let’s define 𝐶 ≡ 𝑎
𝐶 ≡ (𝑎 − 𝑖 𝑏 ) for 𝑛 = 1, 2, 3, … …

𝐶 ≡ (𝑎 + 𝑖𝑏 ) for 𝑛 = 1, 2, 3, … … (3a-c)
where, 𝑖 = −1

By this definition, 𝐶 = 𝐶 (the conjugate of 𝐶 ) (4)

Eq. (3) may be transform into

𝑎 = 𝐶 +𝐶

𝑏 = (𝐶 − 𝐶 ) 𝑖 (5a,b)
Consider an exponential function 𝑒

By Taylor’s expansion,
𝑒 .𝑥 𝑒 .𝑥
𝑒 = 𝑒 + 𝑒 .𝑥 + + +⋯
2! 3!
𝑥 𝑥 (6)
𝑒 = 1 + 𝑥 + + +⋯
2! 3!
Hence,
𝑖𝑥 𝑖𝑥
𝑒 = 1 + 𝑖𝑥 + + +⋯
2! 3!

𝑥 𝑥 𝑥 𝑥
𝑒 = 1 − + −⋯ +𝑖 𝑥 − + −⋯
2! 4! 3! 5!
The Taylor’s expansion of 𝑐𝑜𝑠 𝑥 The Taylor’s expansion of 𝑠𝑖𝑛 𝑥

Euler’s formula 𝑒 = cos 𝑥 + 𝑖 sin 𝑥 (7)

Based on the Euler’s formula, 𝑒 +𝑒
it can be shown that cos 𝑥 =
2

𝑒 −𝑒
sin 𝑥 = (8a,b)
2𝑖

2𝑛𝜋 𝑒 +𝑒 (𝑒 + 𝑒 )
and, in a similar manner, cos 𝑡 = = in a short form
𝑇 2 2

2𝑛𝜋 𝑒 −𝑒 (𝑒 − 𝑒 )
sin 𝑡 = = (9a,b)
𝑇 2𝑖 2𝑖
Introducing equation (5) and (9) into Eq. (1), we obtain

𝑓𝑡 = 𝐶 + 𝐶 +𝐶 𝑒 +𝑒 /2 + 𝐶 −𝐶 𝑒 −𝑒 /2

= 𝐶 + 𝐶 𝑒 + 𝐶 𝑒

𝑓𝑡 = 𝐶 𝑒 (10)

Where
/
1 1 2𝑛𝜋 2𝑛𝜋
𝐶 = 𝑎 −𝑖𝑏 = 𝑓 𝑡 . cos 𝑡 − 𝑖. 𝑓 𝑡 . sin( 𝑡) 𝑑𝑡
2 𝑇 𝑇 𝑇
/
 /
1 1 2𝑛𝜋 2𝑛𝜋
𝐶 = 𝑎 −𝑖𝑏 = 𝑓 𝑡 . cos 𝑡 − 𝑖. 𝑓 𝑡 . sin( 𝑡) 𝑑𝑡
2 𝑇 𝑇 𝑇
/
or
/
𝐶 in the 𝑛𝑡ℎ complex 1
𝐶 = 𝑓 𝑡 .𝑒 𝑑𝑡 (11)
Fourier coefficient 𝑇
/

This equation is valid for 𝑛 = 1, 2, 3, … and 𝑛 = 0 and 𝑛 = −1, −2, −3, …

In polar form:
𝐶 = 𝐶 𝑒 𝐶 is represented by a
point in complex plane.
Fourier amplitude 𝑏
Fourier phase

𝑎
 Amplitude of the 𝑛 harmonic

A periodic function in a discrete function in

time domain frequency domain
Fourier Integral A Fourier transformation for a non-periodic function

Basic Concept: A non−periodic function is a periodic function

where T approaches ∞.

2𝜋
Define ∆𝜔 ≡ (13)
𝑇

Note that when 𝑇 is large, ∆𝜔 is small,

and when 𝑇 approaches ∞, ∆𝜔 becomes infinitely small.

2𝑛𝜋
Define 𝜔 ≡ = 𝑛 ∆𝜔 (14)
𝑇
Let 𝐶(𝜔) can be a continuous function where
/

𝐶𝜔 = 𝑇. 𝐶 = 𝑓 (𝑡) 𝑒 𝑑𝑡 (15)
/
By these definitions, Eq. (10) can be transformed into

∆𝜔
𝑓𝑡 = 𝐶𝜔 𝑒 (16)
2𝜋
Apply the limit,
1 1
lim 𝑓 (𝑡) = lim 𝐶𝜔 𝑒 ∆𝜔 = 𝐶 (𝜔) 𝑒 𝑑𝜔
→ → 2𝜋 2𝜋

and 𝐶 (𝜔) is a continuous function of 𝜔

lim 𝐶(𝜔) = 𝑓 (𝑡) 𝑒 𝑑𝑡
→
In summary, we obtain 1
𝑓𝑡 = 𝐶𝜔 𝑒 𝑑𝜔 (17)
a Fourier Transform pair: 2𝜋

𝐶𝜔 = 𝑓 (𝑡) 𝑒 𝑑𝑡 (18)

in which
𝑓 (𝑡) is a function in time domain, and
𝐶 (𝜔) is a function in frequency domain
𝐶 (𝜔) is called “Fourier amplitude”.

𝐶 (𝜔) is in general a complex function.

Note that one basic requirement for 𝑓 (𝑡) is |𝑓 (𝑡)| 𝑑𝑡 must be finite (not ∞).
 Power Spectral Density Function
1
The Fourier transform of 𝑅 (𝜏) is 𝑆𝜔 = 𝑅 (𝜏) 𝑒 𝑑𝜏 (19)
2𝜋

and Inversely 𝑅𝜏 = 𝑆(𝜔) 𝑒 𝑑𝜔 (20)

where 𝑆 (𝜔) is called the spectral density of the random process.

Eqs (19-20) are called the Wiener-Khintchine Relation. 𝑺(𝝎) and 𝑹(𝝉) are the Fourier transform pair.

Based on some certain properties of 𝑅 (𝜏), it can be shown that 𝑆 (𝜔) is a non-negative function of 𝜔

and 𝑆 −𝜔 = 𝑆(𝜔)
The most important property of 𝑆(𝜔) becomes apparent when we put 𝜏 = 0 in Eq (20):

𝑅𝜏 = 𝑆(𝜔) 𝑒 𝑑𝜔 (20)

𝑅 𝜏=0 = 𝜎 = 𝑆(𝜔) 𝑑𝜔 (21)

The mean square of the stationary random process is given by

the area under a graph of spectral density 𝑆(𝜔) against 𝜔.

𝑆(𝜔) is sometimes called ‘the mean square spectral density’.

Since 𝑅 (𝜏) is a statistical average function (in time domain)

𝑆(𝜔) is also a statistical average function (in frequency domain)

(and it shows the frequency composition of the random process)
 Total area under the curve = 𝜎 = 𝑆(𝜔) 𝑑𝜔

−𝜔 −𝜔 𝜔 𝜔
These two combined areas represent
the contribution of fluctuation of
frequencies between 𝜔 and 𝜔 to the
total mean square 𝜎 .

Wide-band stationary random process. (a) Spectral density and (b) time history of sample Function.
Narrow band and wide band stationary processes
A random process whose spectral density covers a broad band of frequencies is called
a wide band process. Its time history is then made up of the superposition of the whole
band of frequencies.
“Wind velocity” may be treated as a wide band stationary process.

A white noise is an ideal wide band process. It is a physically unrealizable process because
the mean square value of the process is infinite.

But the concept of “white noise” is very useful in the analysis of random response and a spectrum
may be treated as a white noise if it is wide band noise whose bandwidth extends well past all the
frequencies of interest.
A random process whose spectral density covers only a narrow band of frequencies is
called a narrow band process.

In general the response of a lowly

damped structure to a random wind
force is a narrow band process.

The predominant frequency of the

random response will be 𝜔 ,
which is the natural frequency of the
structure.
Spectral Density of 𝑑𝑋/𝑑𝑡 and 𝑑2𝑋/𝑑𝑡2
𝑑𝑋 𝑋(𝑡 + ℎ) − 𝑋(𝑡)
Suppose that 𝑋 (𝑡) is a stationary random process. 𝑋̇ 𝑡 = = lim
𝑑𝑡 → ℎ
𝑋(𝑡 + ℎ) − 𝑋(𝑡) 𝑋(𝑡 + 𝜏 + ℎ) − 𝑋(𝑡 + 𝜏)
𝑅 ̇ (𝜏) = 𝐸 𝑋̇(𝑡) 𝑋̇(𝑡 + 𝜏) = lim 𝐸
→ ℎ ℎ
1
= lim 𝐸 𝑋 𝑡 + ℎ 𝑋 𝑡 + 𝜏 + ℎ − 𝑋 𝑡 + ℎ 𝑋 𝑡 + 𝜏 − 𝑋 𝑡 𝑋 𝑡 + 𝜏 + ℎ + 𝑋(𝑡) 𝑋(𝑡 + 𝜏)
→ ℎ
1
= lim 𝑅 (𝜏) − 𝑅 𝜏 − ℎ − 𝑅 𝜏 + ℎ + 𝑅 (𝜏)
→ ℎ

1 𝑅 𝜏 + ℎ − 𝑅 (𝜏) 𝑅 (𝜏) − 𝑅 𝜏 − ℎ
= lim − −
→ ℎ ℎ ℎ
1 𝑑𝑅 (𝜏) 𝑑𝑅 𝜏 − ℎ 1
= lim − − = lim − 𝑅 (𝜏) − 𝑅 𝜏 − ℎ where 𝑅 ≡
→ ℎ 𝑑𝜏 𝑑𝜏 → ℎ

𝑅 ̇ 𝜏 = −𝑅 (𝜏) where 𝑅 ≡ (22)

𝑅 ̇ (𝜏) is the auto-correlation function o𝑓 𝑋̇(𝑡).

From Eq (20), 𝑅 (𝜏) = 𝑆 (𝜔) 𝑒 𝑑𝜔

𝑑𝑅
= 𝑖 𝜔 𝑆 (𝜔) 𝑒 𝑑𝜔
𝑑𝜏

𝑑 𝑅
= − 𝜔 𝑆 (𝜔) 𝑒 𝑑𝜔
𝑑𝜏

Applying Eq (22) we obtain, 𝑅 ̇ (𝜏) = 𝜔 𝑆 (𝜔) 𝑒 𝑑𝜔 (23)

But since 𝑅 ̇ (𝜏) can also be written as the inverse transform of the spectral density 𝑆 ̇ (𝜔)
we also have
𝑅 ̇ (𝜏) = 𝑆 ̇ (𝜔) 𝑒 𝑑𝜔 (24)
By Comparing Eq (23) and Eq (24), it is clear that

𝑆 ̇ 𝜔 = 𝜔 𝑆 (𝜔) (25)

In similar manner, it can be shown that

𝑆 ̈ 𝜔 = 𝜔 𝑆 (𝜔) (26)

The above two equations can be understood by considering the mean square of a sine function and
its derivatives:
1 1
𝑋 𝑡 = 𝑎 sin 𝜔𝑡 → mean square of 𝑋 = 𝑋 𝑡 𝑑𝑡 = 𝑎
𝑇 2
1 1
𝑋̇ 𝑡 = 𝜔𝑎 cos 𝜔𝑡 → mean square of 𝑋̇ = 𝑋̇ 𝑡 𝑑𝑡 = 𝜔 𝑎
𝑇 2
1 1
𝑋̈ 𝑡 = −𝜔 𝑎 sin 𝜔𝑡 → mean square of 𝑋̈ = 𝑋̈ 𝑡 𝑑𝑡 = 𝜔 𝑎
𝑇 2
Non-stationary random process
Statistical properties of a non-stationary random process change with 𝑡.

Non-stationary in frequency contents

Non-stationary in the mean square value.

𝑆(𝜔, 𝑡) − power spectral density is a function of both 𝜔 and 𝑡

Earthquake ground motion is non-stationary.

Wind velocity fluctuation is approximately stationary.

Harmonic Response of a SDOF System in Complex Form

Consider the equation of motion of a SDOF system, subjected to a unit-amplitude harmonic force:

𝑚 𝑥̈ + 𝑐 𝑥̇ + 𝑘 𝑥 = cos 𝜔𝑡 (a)

or 𝑚 𝑦̈ + 𝑐 𝑦̇ + 𝑘 𝑦 = sin 𝜔𝑡 (b)

Eq (a) + Eq (b) × −1 gives:

𝑚 𝑥̈ + 𝑖 𝑦̈ + 𝑐 𝑥̇ + 𝑖 𝑦̇ + 𝑘 𝑥 + 𝑖 𝑦 = cos 𝜔𝑡 + 𝑖 sin 𝜔𝑡

𝑚 𝑧̈ + 𝑐 𝑧̇ + 𝑘 𝑧 = 𝑒 (c)

a unit amplitude complex excitation

where 𝑧 𝑡 = 𝑥 𝑡 + 𝑖 𝑦(𝑡)

𝑧(𝑡) is called the “complex response”.

In steady state, it can be shown that
𝑚 𝑥̈ + 𝑐 𝑥̇ + 𝑘 𝑥 = cos 𝜔𝑡
𝑥 𝑡 = 𝑎 cos 𝜔𝑡 − 𝜙
𝑚 𝑦̈ + 𝑐 𝑦̇ + 𝑘 𝑦 = sin 𝜔𝑡
𝑦 𝑡 = 𝑎 sin 𝜔𝑡 − 𝜙 (d)

( )
𝑧 𝑡 = 𝑥 𝑡 + 𝑖 𝑦 𝑡 = 𝑎 cos 𝜔𝑡 − 𝜙 + 𝑖. 𝑎 sin 𝜔𝑡 − 𝜙 = 𝑎 𝑒 = 𝑎𝑒 .𝑒
(e)
𝑎𝑒 is called “complex amplitude” and it is a function of 𝝎.

We may write Eq (e) in a general form as: 𝑧 𝑡 = 𝐻 𝜔 .𝑒 (f)

where 𝐻 (𝜔) is the complex amplitude or the “complex frequency response function”.

𝑅𝑒𝑠𝑝𝑜𝑛𝑠𝑒 𝑂𝑢𝑡𝑝𝑢𝑡
𝐻 𝜔 = =
𝐹𝑜𝑟𝑐𝑒 𝐼𝑛𝑝𝑢𝑡

We can construct 𝐻 (𝜔) from 𝑎 (𝜔) and 𝜙 𝜔 since 𝐻 𝜔 = 𝑎 𝑒

Alternatively, we may substitute Eq (f) into Eq (c): 𝑧 𝑡 = 𝐻 𝜔 . 𝑒 𝑚 𝑧̈ + 𝑐 𝑧̇ + 𝑘 𝑧 = 𝑒

𝐻 (𝜔) −𝑚 𝜔 + 𝑖𝑐𝜔 + 𝑘 𝑒 =𝑒

1
That is 𝐻𝜔 = (g)
𝑘 − 𝑚𝜔 + 𝑖𝑐𝜔
1
Therefore 𝑧 𝑡 = 𝐻 𝜔 .𝑒 = 𝑒
𝑘−𝑚𝜔 +𝑖𝑐𝜔

1
or 𝑧𝑡 = 𝑒 (h)
𝑚 𝜔 −𝜔 +2𝜉𝜔𝜔 𝑖
𝐻 (𝜔)

where 𝜔 = 𝑘/𝑚 and 𝜉 = 𝑐/(2 𝑘 𝑚 )

Note that 𝑎(𝜔) = 𝐻(𝜔) and 𝜙(𝜔) = −Arg 𝐻 𝜔 (i)

Response of a SDOF system to a stationary random force

Equation of motion 𝑚 𝑥̈ + 𝑐 𝑥̇ + 𝑘 𝑥 = 𝑓(𝑡) (a)

Response: 𝑥 𝑡 = 𝑥 𝑡 + 𝑥 (𝑡)

The homogenous solution: 𝑥 𝑡 →

0 in the stationary state of response

Hence, in the stationary state 𝑥𝑡 = 𝑥 𝑡 = 𝑓 (𝜏) ℎ(𝑡 − 𝜏) 𝑑𝜏 (b)

where ℎ(𝑡) is the unit impulse response function and is given by

1
𝑒 sin 𝜔 𝑡 for 𝑡 > 0
ℎ(𝑡) = 𝑚𝜔 (c)
0 for 𝑡 ≤ 0
Since ℎ 𝑡 − 𝜏 = 0 for 𝜏 > 𝑡 , equation (b) can be transferred into

𝑥𝑡 = 𝑓 (𝜏) ℎ(𝑡 − 𝜏) 𝑑𝜏

By changing the variable 𝜏 to 𝜃 where 𝜃 ≡ 𝑡 − 𝜏 , we obtain

𝑥𝑡 = 𝑓 (𝑡 − 𝜃) ℎ(𝜃) 𝑑𝜃 (d)

For a special case where 𝑓 𝑡 = 𝑒

𝑥 (𝑡) = 𝑒 ( ) ℎ(𝜃) 𝑑𝜃
For a special case where 𝑓 𝑡 = 𝑒

( )
𝑥𝑡 = ∫ 𝑒 ℎ(𝜃) 𝑑𝜃

𝑥𝑡 = ℎ𝜃 𝑒 𝑑𝜃 . 𝑒

From the previous section, we got 𝑥(𝑡) = 𝐻 (𝜔). 𝑒

By comparing the above two equations, it is clear that

𝐻𝜔 = ℎ𝜃 𝑒 𝑑𝜃 (e)

Equation (e) says that the complex frequency response function 𝐻 (𝜔) is the Fourier transform
of the unit impulse response ℎ(𝑡).
Next let’s consider the case where 𝑓 (𝑡) is a stationary random force.

𝑋(𝑡) = 𝑓 (𝑡 − 𝜃) ℎ(𝜃) 𝑑𝜃
A stationary random response A stationary random force

The mean response to random force is 𝜇 𝑡 = 𝐸 𝑋 𝑡 =𝐸 𝑓 (𝑡 − 𝜃) ℎ(𝜃) 𝑑𝜃 (f)

It follows that 𝜇 (𝑡) = 𝐸 𝑓 (𝑡 − 𝜃) ℎ(𝜃) 𝑑𝜃

S 𝐸 𝑋+𝑌 =𝐸 𝑋 +𝐸 𝑌
𝜇 (𝑡 − 𝜃) = 𝜇

𝐸 𝐶𝑋 = 𝐶 𝐸𝑋
Hence 𝜇 = 𝜇 ℎ(𝜃) 𝑑𝜃 (g)
 𝜇 = 𝜇 ℎ(𝜃) 𝑑𝜃 (g)

1
It is clear from Eq (e) and Eq (g) that 𝐻 𝜔 = ℎ 𝜃 𝑒 𝑑𝜃 𝐻 𝜔 =
𝑘 − 𝑚𝜔 + 𝑖𝑐𝜔

1
ℎ(𝜃) 𝑑𝜃 = 𝐻 𝜔=0 = (h)
𝑘

𝜇
Hence 𝜇 = (i)
𝑘

For the case of zero-mean random force, 𝜇 =0 → 𝜇 =0

Zero mean force → zero mean response
The autocorrelation function of the random response, denoted by 𝑅 (𝜏) is then given by

𝑅 𝜏 = 𝐸 𝑋 𝑡 𝑋 𝑡+𝜏 = 𝐸 𝑓 (𝑡 − 𝜃 ) ℎ(𝜃 ) 𝑑𝜃 𝑓 𝑡 + 𝜏 − 𝜃 ℎ 𝜃 𝑑𝜃

= 𝐸 ℎ 𝜃 ℎ 𝜃 𝑓(𝑡 − 𝜃 ) 𝑓 𝑡 + 𝜏 − 𝜃 𝑑𝜃 𝑑𝜃

= ℎ 𝜃 ℎ 𝜃 𝐸 [𝑓 (𝑡 − 𝜃 ) 𝑓 𝑡 + 𝜏 − 𝜃 ] 𝑑𝜃 𝑑𝜃

𝑅 𝜏 = ℎ 𝜃 ℎ 𝜃 𝑅 𝜏−𝜃 +𝜃 𝑑𝜃 𝑑𝜃 (j)

Double integration—not convenient to use for analysis purpose

The spectral density function of the random response, denoted by 𝑆 (𝜔) is given by

1
𝑆 𝜔 = 𝑅 𝜏 𝑒 𝑑𝜏 (k)
2𝜋

Eq (k) is the Weiner-Khinchin relation, as described earlier by equation (19).

Introducing Eq (j) into Eq (k), we obtain:

1
𝑆 𝜔 = ℎ 𝜃 ℎ 𝜃 𝑅 𝜏−𝜃 +𝜃 𝑒 𝑑𝜃 𝑑𝜃 𝑑𝜏
2𝜋

1
= ℎ𝜃 𝑒 ℎ𝜃 𝑒 𝑅 𝜏−𝜃 +𝜃 𝑒 𝑑𝜃 𝑑𝜃 𝑑𝜏
2𝜋

By changing the variable 𝜏 to 𝜃 where 𝜃 ≡ 𝜏 − 𝜃 + 𝜃 , it follows that

 1
𝑆 𝜔 = ℎ𝜃 𝑒 𝑑𝜃 ℎ𝜃 𝑒 𝑑𝜃 𝑅 𝜃 𝑒 𝑑𝜃
2𝜋

= 𝐻 −𝜔 . 𝐻 (𝜔). 𝑆 (𝜔) see Eq (e)

Since 𝐻 (−𝜔) is the complex conjugate of 𝐻 (𝜔), 𝐻 −𝜔 . 𝐻 𝜔 = 𝐻 (𝜔)

Therefore, 𝑆 𝜔 = |𝐻 (𝜔)| . 𝑆 (𝜔) (l)

Simple algebraic relation

1
Where 𝐻𝜔 =
𝑚 𝜔 −𝜔 +2𝜉𝜔𝜔 𝑖
1
and 𝐻 (𝜔) = (m)
𝑚 𝜔 −𝜔 +4𝜉 𝜔 𝜔
 1
𝐻(𝜔) 𝐻 (𝜔) =
𝑚 𝜔 −𝜔 +4𝜉 𝜔 𝜔

1
4𝜉 𝑘

1 1
=
𝑘 𝑚 𝜔

𝜔 𝜔
 𝑓 (𝑡): a stationary random force

This provides a physical insight into the

characteristics of random excitations
and responses
A SDOF system (which acts like a mechanical filter)

𝑋(𝑡): a stationary random response

Mean Square of Response

Root-mean- square of response: 𝜎 = 𝑆 (𝜔) 𝑑𝜔 = 𝐻𝜔 . 𝑆 𝜔 𝑑𝜔 (a)

For a special case where 𝑓 (𝑡) is a stationary white noise with spectral intensity So, that is,
𝑆 𝜔 = 𝑆

𝜎 = 𝑆 𝐻𝜔 𝑑𝜔 (b)

The integral in this R.H.S can be evaluated by means of the complex integral using the residual
at the pole, that is,
𝑖𝜔𝐵 +𝐵 where 𝐴 , 𝐴 , 𝐴 , 𝐵 , 𝐵
for 𝐻𝜔 =
−𝜔 𝐴 + 𝑖 𝜔 𝐴 + 𝐴 are real constants

𝐵 /𝐴 𝐴 + 𝐵
|𝐻 (𝜔)| 𝑑𝜔 = 𝜋 (c)
𝐴 𝐴
 𝜋𝑆 𝜋 𝑆 𝜔 1
Applying Eq (c), 𝜎 = = (d)
2𝑚 𝜔 𝜉 2 𝑘 𝜉

For a realistic case where 𝑆 (𝜔) is rather uniform but not “white” and the structural system
is lightly damped:
𝜋 𝑆 (𝜔 )
𝜎 ≈ (e)
2𝑚 𝜔 𝜉
Because the contribution of PSD around 𝜔 is the most significant.

Note that the rms of random response 𝜎 ∝ 1/ 𝜉

While the response amplitude to a sine function (at resonant peak) ∝ 1/ 𝜉

Extremes in Narrow-band Random Response

Let’s consider two basic types of structural failure:

1) Failure due to first response excursion up to structural limit.

2) Failure due to an accumulation of damage.
 Peak: an extreme of 𝑋 (𝑡): a candidate for maximum response.

Assuming that the response is a stationary narrow-band Gaussian

process, it can be shown that

∫ 𝜔 𝑆 𝜔 𝑑𝜔
(a) Statistical average frequency “𝜔 ” 𝜔𝒂𝟐 =
∫ 𝑆 𝜔 𝑑𝜔
𝜔 is simply a weighted average of 𝜔
in which 𝑆(𝜔) is the weighting function.

It can be shown that the frequency of “zero-crossing

with positive slope”, 𝜆, is equal to 𝜔 /2𝜋.

For a narrow band response of a SDOF structure with natural frequency 𝜔 , 𝜔 ≈ 𝜔

The occurrence rate of peak 𝜈 ≈ 𝜆 ≈ 𝜔 /2𝜋

(b) Probability density function of peak 𝑓𝑝(𝑎): Peak amplitude (p) is a random variable !
𝑋(𝑡) 𝑋(𝑡)
a narrow band random process a wide band random process

𝑡 𝑡

𝑓𝑝(𝑎) 𝑓𝑝(𝑎)
𝑎 1
𝑒 𝑒
𝜎
2𝜋𝜎

𝑎 𝑎
𝜋
𝜇 = 𝜎 𝜇 =0
2
Rayleigh distribution Gaussian distribution
In reality, 𝑓𝑝(𝑎) lies between these two extremes.

𝑓𝑝(𝑎)
Real condition of PDF
(depends on spectral shape)

𝑎
(c) Distribution of maximum peak:
Consider a response in time duration 𝑇, let’s define 𝑎 = 𝑚𝑎𝑥 {𝑎 , 𝑎 , 𝑎 … . 𝑎 }

𝑓𝑝(𝑎)

Let’s assume that 𝑓𝑝(𝑎) is known.

Now, we want to know the probability density function of 𝑎

𝑛: the number of peaks in the duration 𝑇

𝐸 𝑛 = 𝜈𝑇

Using “the statistical theory of extremes”,

it can be shown that.

𝜇 = 𝑔𝜎

Where 𝑔 is called ‘peak factor’ and is given by

0.577
𝑔 = 2 ln(𝜈𝑇) +
2 ln 𝜈𝑇

𝜇 is the mean value of “maximum” peak within time duration 𝑇

Numerical Example:

𝜈 = 1 peak/ sec (𝜔 /2𝜋 = 1 Hz)

𝑇 = 1000 sec

𝜈𝑇 = 1000

𝜇 = average of 𝑎 = 3.9 𝜎

The maximum peak within 1000 sec would be

about four times of the root-mean- square of
the response.
 d) Fatigue life
S: Stress amplitude
Conventional fatigue test
N: Number of cycles to failure

Results are typically c and b are constants where c, b > 0

presented in the form of
SN curves These values depend on the type of
material, welding characteristics, level
of mean stress etc.
Problem:

Conventional test condition: constant amplitude stress S

Real condition: random amplitude stress S How to predict fatigue life? (based on the test data)
Palmgram and Miner’s “Linear Damage Accumulation Hypothesis”

When 𝑛 cycles of stress amplitude 𝑆 have been experienced,

the material has “used up” a fraction of its fatigue life equal to
𝑛/ 𝑁 where 𝑁 is the number of cycles at which failure occurs
under uniform amplitude 𝑆 as indicated by SN curve.

Fractional damage 𝑑𝑖 due to stress amplitude 𝑆 and 𝑛 cycles:

𝑛
𝑑 = where 𝑁 =𝑐𝑆
𝑁

𝑛
Total damage: 𝐷𝑇 = 𝑑𝑖 =
𝑁

Fatigue failure occurs when 𝐷𝑇 > 1

For a narrow-band random response:
𝑓 𝑎
𝑎 𝑎
𝑓 𝑎 = exp −
𝜎 2𝜎

Time duration of response = 𝑇

Total number of peaks in 𝑇 = 𝜈𝑇 Peak distribution of stress: Rayleigh distribution

(𝑎 is stress peak or stress amplitude)

The expected number of peaks which lie

𝑛 𝑎 = 𝜈𝑇. 𝑓 (𝑎) 𝑑𝑎
between 𝑎 and 𝑎 + 𝑑𝑎 = 𝑛(𝑎)

The expected fractional damage: 𝑛(𝑎) 𝑓 (𝑎) 𝑑𝑎

= 𝜈𝑇 where 𝑁 𝑎 =𝑐𝑎
𝑁(𝑎) 𝑁(𝑎)
𝑓 (𝑎)
The total expected damage 𝐷𝑇 ∶ 𝐷 = 𝜈𝑇 𝑑𝑎
𝑁(𝑎)
 𝑓 (𝑎) 𝑓 𝑎
𝐷 = 𝜈𝑇 𝑑𝑎
𝑁(𝑎)

𝜈𝑇
𝐷 ≅ 2𝜎 .Γ 1 + 𝑏 2
𝑐

where Γ : gamma function Peak distribution of stress: Rayleigh distribution

(𝑎 is stress peak or stress amplitude)

Γ 𝜃 =2 𝑥 exp −𝑥 𝑑𝑥

The expected damage accumulation increases in linear

proportion to the number of cycles 𝜈𝑇. It also depends
(nonlinearly) on the rms level (𝜎 ) of the stress history
and on the constants 𝑏 and 𝑐.

When 𝐷 =1 fatigue failure

 𝜈𝑇
𝐷 ≅ 2𝜎 .Γ 1 + 𝑏 2
𝑐

By setting 𝐷 = 1, T becomes 𝑇 ∗

𝑐 1
𝑇∗ = 2𝜎 .Γ 1+𝑏 2
Γ (𝑥) =
𝜈 Γ(𝑥)

𝑇 ∗ : average fatigue life

𝑐, 𝑏: constants obtained from material testing

𝜈, 𝜎 : response parameters obtained from

random vibration analysis
Numerical example:

75S-T6 aluminum alloy

𝑏 = 6.09 obtained from material testing,
. where stress is expressed in “psi” unit.
𝑐 = (2 × 10 ) = 1.92 × 10

3000
𝜎 = 𝑝𝑠𝑖
𝜉 obtained from random vibration analysis.
𝜈 = 1 cycle/ sec (Hz) Γ 𝑥

6.09
Γ 1+ = 6.35
2
.
∗
1.92 × 10 2 . 3000 1
𝑇 = × ×
1 𝜉 6.35

𝑇 ∗ = 2.44 × 10 (𝜉) .
𝑥
 Fatigue Life

𝑇 ∗ = 2.44 × 10 (𝜉) .

𝜉 𝜎 (𝑝𝑠𝑖) 𝑇 ∗ (𝑠𝑒𝑐) 𝑇 ∗ (ℎ𝑟)

0.01 30000 1938 0.54

0.02 21213 16052 4.5

0.05 13416 262572 73

Thank You

77