Department of Mechanical Engineering Blekinge Institute of Technology

Master's Degree Thesis ISRN: BTH-AMT-EX--2007/D-07--SE

Supervisor: Kjell Ahlin, Professor Mech. Eng.

Department of Mechanical Engineering Blekinge Institute of Technology

Karlskrona, Sweden 2007

Ejaz Yousaf

Output only Modal Analysis

Test Structure

?

Scaled

Mode Shapes

Output only modal Analysis

Ejaz Yousaf

Department of Mechanical Engineering Blekinge Institute of Technology

Karlskrona, Sweden 2007

Thesis submitted in fulfillment of the requirements for the award of the Master of Science in Mechanical Engineering by the Department of Mechanical Engineering, Blekinge Institute of Technology, Karlskrona, Sweden.

ABSTRACT

This thesis work aims at estimating the modal parameters of a system using a technique amongst different proposed in Output only modal analysis. This technique assumes that impulse response function can be approximated by cross correlation between responses and transfer function by cross power spectral density of the responses. A method which involves added masses on some points of structure under test is used to obtain a scaling factor. This scaling factor is used for scaling the unscaled modal vectors from Output- Only modal testing. A simulation on four degrees of freedom system in Matlab is performed in order to check the validation of theory. Experiments are performed on a structure with bolt channel joints. Modal identification tool box developed by Saven Edutech AB in Matlab and a few Matlab built-in functions are used during the analysis of the thesis work. All the structures under study in this thesis are modeled as linear time invariant systems.

Keywords

Cross correlation function, Cross power spectral density, FRF, Impulse response

function, Scaling factor, Scaled mode shapes, ambient excitations

Acknowledgements

Some people are more respectful and nice than what they look to us. Their personality seems far better than we describe in our words. There are a number of personalities, who have helped and contributed in what I have achieved today and I am grateful to all of them.

Saying the word “Thanks” to someone means paying back for what they have given to you or helped you in achieving something, but I cannot pay back to Professor Kjell Ahlin and my parents just by saying “Thanks” for their guidance, help and supporting me throughout my career especially my Masters studies. I will remember their contributions, guidance and help for the rest of my life. It would have been impossible for me to reach up to this stage without their support and guidance. I am very grateful to all of them.

I would like to mention the names of Mr. Martin Magnevall and Andreas Josefsson (PhD.

Students) for their help, guidance and support throughout the course work and thesis.

I am also grateful to the whole department of Mechanical Engineering at Blekinge Institute of Technology, Karlskrona and to all those, whose published data and materials are used as references in this thesis work.

Karlskrona April, 2007

Ejaz Yousaf

Contents

1 NOTATIONS ... 5

1.1 L IST OF SYMBOLS ... 5

1.2 L IST OF OPERATORS ... 6

1.3 L IST OF ABBREVIATIONS ... 7

2 INTRODUCTION ... 8

2.1 B ACKGROUND ... 8

2.2 E XPERIMENTAL M ODAL A NALYSIS ... 8

2.3 S TRUCTURAL M ODIFICATION ...10

2.4 O UTPUT O NLY M ODAL A NALYSIS ...11

2.4.1 N EED ...11

2.4.2 A PPLICATIONS ...11

2.4.3 P ROCEDURE ...12

2.4.4 T HINGS TO R EMEMBER ...15

3 OUTPUT ONLY SYSTEM IDENTIFICATION ...16

3.1 S YSTEMS I DENTIFICATION OF C IVIL E NGINEERING S TRUCTURES ...16

3.1.1 N ON P ARAMETRIC METHODS ...17

3.1.2 P ARAMETRIC METHODS ...17

4 PARAMETER IDENTIFICATION METHOD ...20

4.1 P OWER S PECTRAL D ENSITY E STIMATION ...20

4.1.1 P ERIODOGRAM M ETHOD ...20

4.1.2 C ORRELATION M ETHOD ...22

4.2 M ODAL D ECOMPOSITION OF S PECTRA D ENSITIES ...22

4.3 T RUNCATION OF THE C ROSS -S PECTRAL D ENSITY M ATRIX ...25

4.4 P OLES AND R ESIDUES E STIMATION ...25

4.5 S ELECTION OF R EFERENCE R ESPONSES ...26

5 NORMALIZATION OF MODE SHAPES ...27

5.1 DERIVATION ...27

5.1.1 S CALING FACTOR ...30

5.2 T HINGS TO R EMEMBER ...34

6 SIMULATION...35

6.1 D ESCRIPTION OF S YSTEM ...35

7 APPLICATION TO A BOLT CHANNEL JOINT STRUCTURE...41

7.1 I NTRODUCTION ...41

7.2 E QUIPMENT S ETUP ...41

7.3 E XPERIMENTAL S ETUP ...43

7.3.1 P ROCEDURE ...44

7.3.2 E STIMATION OF M ODE S HAPES ...47

7.4 C OMPARISON OF R ESULTS ...48

7.5 D ISCUSSIONS ...53

CLOSING REMARKS ...56

REFERENCES ...57

1 Notations

1.1 List of symbols

F(t) excitation force at time t f Frequency [Hz]

f

s

Sampling frequency [Hz]

(ω )

H Frequency response function I Identity matrix

j Imaginary unit j

²

= − 1 k Discrete time instant t = k t Δ M, C, K Mass, damping and stiffness matrix

N

i

Number of inputs

N

o

Number of outputs

N

m

Number of DOF

N

ref

Number of references N

it

Number of iterations

N

a

Number of averages

Q Modal participation factor

R Residue matrix

R

i

Output covariance matrix at time lags i R

xy

Covariance between variables x and y r Number of references

(ω )

S Power spectrum continuous time variable u(t) Input at time t

u

k

Input at time k

V Mode shape matrix V = C Ψ

x(t) State at time t

y(t) Output at time t

y

k

Output at time instant k

ref

y

k

Reference output at time instant k

ξ Damping ratio

ψ Operational mode shape

φ Normalized mode shape

α Operational scaling factor

Δ m Local change in mass

Δ t Sampling interval

)

δ (t Dirac delta

δ Kronecker

k

delta

λ

i

Continuous-time Eigen value μ

i

Discrete-time Eigen value

ω Angular frequency [rad/s]

1.2 List of operators

)

T

(⋅ Matrix transpose

)

1

( ⋅

⁻

Matrix inverse )

*

(⋅ Complex conjugate

)

H

(⋅ Hermitian transpose: complex conjugate transpose of matrix Re ) (⋅ Real part of

Im ) (⋅ Imaginary part of

1.3 List of abbreviations

AR Autoregressive

ARMA Autoregressive Moving Average CMIF Complex Mode Indicator Function DFT Discrete Fourier Transform DOF Degree of Freedom EMA Experimental Modal Analysis EMT Experimental Modal Testing

FE Finite Element

FEM Finite Element Method FFT Fast Fourier Transform FRF Frequency Response Function IFFT Inverse Fast Fourier Transform

LS Least Squares

LSCE Least-Square Complex Exponential MAC Modal Assurance Criterion ODS Operational Deflection Shape OMA Operational Modal Analysis OMT Output-only Modal Testing

PP Peak Picking

RMS Root Mean Square

SNR Signal to Noise Ratio

SSI-COV Covariance-driven Stochastic Subspace Identification SSI-DATA DATA-driven Stochastic Subspace Identification SVD Singular Value Decomposition

RD Random Decrement

RMITD Reference Ibrahim Time Domain

ITD Ibrahim Time Domain

2 Introduction

2.1 Background

Main target of Engineers has been to make the end product qualitative and cost effective with short lead time. This phenomenon lead to improvement in the technological methods involved in the testing and design phases of a product.

The dynamical systems vibrate during operation as the body does not remain in equilibrium. If the frequency of excitation force matches with any resonance frequency of the structure, it could vibrate with unexpectedly high amplitudes causing damage or failure of the structure.

Modal Analysis is an important tool in vibration analysis of a structure. Objective of modal analysis is the evaluation of modal parameters of a structure under dynamic loads.

Modal parameters provide information of protecting the structure from failure.

Resonance frequencies, damping ratios and mode shapes are the required modal parameters of a dynamical system. Such frequencies at which structure vibrates violently even at small amplitude excitation are known as resonance frequencies. Resonance frequencies are the Eigen values of the equation of motion. Damping ratios decide how quickly or slowly the structure can dissipate the vibration energy. The relative position of points on a structure at a given natural frequency is known as mode shape. Mode shapes are termed as Eigen vectors.

2.2 Experimental Modal Analysis

Finite element analysis FEA uses the Eigen value solvers to estimate the modal parameters. When experiment is performed on real life structure, we divide the structure into discrete points. But true information regarding mass, stiffness and damping matrices is missing for the chosen points. In Experimental modal analysis, we excite the structure (location depend on the selection of excitation) and measure the response at all the points.

Different methods are available in both time and frequency domain to evaluate the modal parameters. Measurement data is usually used to calculate the impulse response and frequency response functions. Modal parameters are extracted from these two functions.

A frequency response function is estimated by the equation (2.1)

) (

) ) (

( G f

f f G

H

XX

=

YX

(2.1)

Where, ) ( f

G

_XX

is the auto power spectral density of the input signal )

( f

G

_YX

is the cross power spectral density of the input and output signal )

( f

H is the frequency response function

Steps used in traditional experimental modal analysis are presented in the following block diagram.

Figure 2.1 Block Diagram of Experimental Modal Analysis

2.3 Structural Modification

After agreement between experimental modal analysis and finite element analysis, our modal model is developed. Model is simulated in the operating conditions. If model operates safely under the operating conditions, it marks the end of the modal analysis design phase. If model is not operating safely in the operating conditions, the structural dynamic modification (SDM) is carried out. Modification is done until required conditions are satisfied.

Figure 2.2 Design Phases of Modal Analysis and Structural Dynamic Modification

Experimental Modal Testing

Finite Element Modeling

Modal Parameter Estimation

Eigen Solution

Development of Modal Model

Structural Changes required?

Use SDM

Yes

DONE

No

2.4 Output Only Modal Analysis

In traditional modal analysis, excitations and responses of the structure are measured to extract the modal parameters. If designer could extract the modal parameters with as high accuracy as with traditional modal analysis without measuring the input force, it will decrease the project cost as well as analysis time. Output only modal analysis is a technique in which we do not need to measure the input excitations and just responses are measured. Assumption in output only modal analysis is that the excitation force is white noise. Output-Only modal analysis is also known as Operational Modal analysis because it deals with testing of structure during normal operation.

2.4.1 Need

For large engineering structures like bridges, dams and high rise buildings, it is difficult to excite them artificially and to measure the excitations. Random forces act on them such as on bridge different forces such as vehicle trafficking, wind excitation, ocean waves etc excite the structure together so it’s almost impossible to measure all these forces simultaneously. If the forces are not measured correctly, then modal analysis cannot give accurate estimates of the modal parameters. On the other hand, when we excite such structure, large amplitude force is required to vibrate the structure at all points under observation. Such force can cause local damages to the structure. Another solution to the problem was advised, which state that model of structure should be tested in laboratory to estimate the modal parameters. Natural conditions under which structure operates are very difficult to be produced in the laboratory. So, output only modal analysis is the best option available in these conditions as it depends on the responses only and the responses could be measured with high accuracy.

Output only modal analysis is good in a way that we do not need to make laboratory models for it. We can perform the experiments on the site of the structure. It means that structure is not disturbed of its usage and modal analysis could be performed during the normal operating conditions of the structure. To obtain a very accurate modal model of a structure, huge data is required which means long testing time is required i.e. for days or weeks. Traditional modal analysis needs to stop the operations of the structure in case of monitoring. On the other hand, output only modal analysis considers the natural excitations as input so all we need are the installation of sensors.

2.4.2 Applications

Output only modal testing is used in many engineering structures such as automotives,

bridges, high rise buildings, space shuttles, dams etc.

Like traditional modal analysis, finite element model updating is done with output only modal analysis. Updated finite element model can predict the dynamic behavior of the structure more reliably.

Obtained modal model of different components of complex structures the dynamical behavior of the complete structure can be computed by using sub-structuring techniques [1].

A structure suffering from vibration problem could be modified by structural dynamics modification. So modal model obtained can be modified by making changes in mass, stiffness and damping at suitable locations in order to get rid of unwanted vibrations.

In the stage of product design and development, when designers agree with a model, it is considered as a reference modal with known modal parameters. During operation, structural health monitoring is done in order to check faults, performance degradation, component deterioration, impending failure etc. [2]. So, the standard reference model is compared with in-operation model to predict the damages. Non destructive testing and fatigue analysis are the concerned fields in structural health monitoring [3].

If the vibration level needs to be documented in locations where no measurements can be made, Output only modal analysis can do it if we have a Finite Element Model available.

From the modal test we will obtain the modal coordinates or modal response at some measurable locations. These modal responses will then be extrapolated to other unmeasured locations through the mode shapes of the FE model, and by superposition the actual responses at the location is estimated. Even though the FE models only return normal modes this extrapolation will be good most structures.

2.4.3 Procedure

In Output-Only modal analysis, there are different techniques followed to estimate the modal parameters. A discussion on different methods is detailed in chapter 2 followed by the theoretical treatment and aspects of the techniques in chapter 3.

Method used in the thesis work approximates FRFs by the cross PSD between the responses and Impulse response functions by the cross correlations between responses.

Details of the relevant theory are discussed in chapter 3 and 4.

To calculate the modal participation factor, force exciting the structure should be

measured. This is the requirement to obtain scaled mode shapes. In output only modal

analysis, only responses are measured which means that scaled mode shapes could not be

obtained like in traditional modal analysis. During past few years, many techniques have

been developed to obtain scaled mode shapes from output only modal analysis.

One simple method amongst devised methods is mass changes in different locations of the structure. A scaling factor which is equivalent to the modal participation factor is obtained by this technique formulated in chapter 5. Using this scaling factor, unscaled mode shapes obtained from response data are scaled.

Experiment is held in two stages. In first stage, modal parameters are extracted on the reference structure which is not loaded by masses. Same procedure is repeated on the structure loaded by masses at certain locations. Modal parameters obtained from the two stages experiment are used to derive a scaling factor. Details of derivation of scaling factor are given in Chapter 5. Mode shapes are scaled using this scaling factor.

Outlines of the procedure employed in the thesis work are shown in the block diagram,

figure 2.3. One location is considered as reference. Cross correlations and cross PSD’s

are calculated between the reference and rest of locations. Methods used to estimate the

modal parameters are discussed in 3 and 4.

Figure 2.3 Procedure of Estimation of Modal Parameters employed in thesis work

Test Structure

Response from Reference location

Response from Rest of locations

Cross Correlations b/w responses

Cross PSDs b/w responses

MODAL PARAMETERS

• Resonance Frequencies

• Damping

• Modal Vector (unscaled)

Repeat experiment with added masses on structure

MODAL PARAMETERS

• Resonance Frequencies

• Damping

• Modal Vector (unscaled)

Scaling Factor

Scaled

Mode Shapes

Initially, theory is applied in a simulation on a four DOF system detailed in chapter 6.

Proposed method is applied by performing required experiments to a bolt channel joint structure detailed in chapter 7.

2.4.4 Things to Remember

In Output only modal testing, from experimental view we have to take care of a few general things.

1 Structure is divided in to discrete locations where responses are to be measured.

There are certain locations on a structure where structure does not vibrate. Such locations are known as nodes. These locations should be avoided. In order to be sure that enough and correct locations are selected, a finite element model is useful in the initial stages.

2 We should get good data. Signal to noise ratio must be good and we should remove spikes from data. Some time domain techniques are applied to the collected data which give smoothness to the data.

3 We should collect enough data in order to get nice frequency resolution and to get more averaging with large block sizes. Bias errors are removed by selecting large block sizes, while increasing the averages reduce the random error in the spectrums.

4 Excitation signal should have a flat spectrum in the vicinity of resonances of the

structure. Unevenness in the excitation spectrum produces skewness in the frequency

response function.

3 Output Only System Identification

Mathematical models are a convenient way of describing a dynamic system.

Mathematical models are described in two ways [4]

● Physical modeling

● System Identification

In physical modeling, dynamic model is constructed on physical knowledge and fundamental laws of motion i.e. M, C and K matrices should be known. If physical knowledge is limited, a model based on input/output behavior of the system can be obtained through system identification based on calibration of a model using experimental data. System identification is not a complete substitute of physical modeling. System Identification basically means modeling of a dynamic system from experimental data. Our aim in system identification is to produce modal model.

3.1 Systems Identification of Civil Engineering Structures

System is excited by Impulse hammers or electro dynamic shakers while responses are recorded usually by accelerometers. According to traditional modal analysis, inputs to the system i.e. excitation forces must be measured to estimate a perfect modal model.

Measurement of excitation is very easy when the excitation force is the dominant and noise is negligible. In most of the cases, noise is of negligible level as compared to the input excitation. In such cases, we can predict modal model which is reliable as changes in the modal model due to disturbance are negligible.

Excitations on large structures are difficult to measure. Also, exciting these structures with artificial exciters is a tough job because of many practical difficulties. Like on a bridge, wind excitation, vehicle trafficking, water waves are the three primary excitations. It’s difficult to measure all of them as the nature of excitation is random at different locations. During such a case, we cannot measure the excitations accurately and it results in an erroneous modal model.

The only parameter which could be measured accurately is the response data. Output

only modal analysis is system identification method which estimates a modal model by

using responses data only. Input excitation is considered to be white noise.

There are two main groups of output-only modal identification methods: nonparametric methods essentially developed in frequency domain and parametric methods in time domain .

3.1.1 Non Parametric methods

Following non parametric methods are used in practice.

3.1.1.1 Peak Picking Method

In this method, average normalized power spectral densities (ANPSDs) and ambient response functions between all the measurement points of the structure are evaluated.

Peaks in the spectra give the Eigen frequencies and modal damping is obtained by half power band width method. For more details about half power band width, please see [5].

This method has been used in the modal identification of buildings [6], [7] and bridges [8], [9]. Felber [10] systematized this method later.

3.1.1.2 Frequency Domain Decomposition

Power spectral densities of a set of SDOF systems are evaluated by performing single value decomposition of the response spectra matrix. This method was systematized and detailed by Brincker [11]. In extracting modal damping factors estimate, more enhancements were done by [12].

3.1.1.3 Enhanced Frequency Domain Decomposition

In this method, first we obtain the power spectral densities of a set of SDOF systems.

Then the auto correlation functions are evaluated by taking inverse Fourier transform of power spectral densities. Inspections of modal parameters are done from the inspection of the decay of auto correlation functions.

3.1.2 Parametric methods

Parametric methods for modal identification are as follows.

3.1.2.1 Stochastic Subspace and Auto regressive moving averages (ARMA) models

An idealized mathematical model for the structural behavior of the system is devised in

these methods. These methods can be directly applied to response time data or to

response correlation functions. These functions are either evaluated using FFT algorithm

[13] or by using Random Decrement method (RD) [14]. Modal parameters are estimated

and the model is fitted to experimental data as much as possible. For more details about

Stochastic subspace method and ARMA models, interested reader is referred to [15], [16]

and [4].

3.1.2.2 Least Squares Complex Exponential, Poly Reference Complex Exponential and Covariance- Driven Stochastic Subspace Identification methods These methods use response correlation functions to estimate the modal parameters. See [17] for Least squares complex exponential (LSCE), [18] for Poly reference complex exponential (PRCE), [15] for covariance driven stochastic subspace identification methods for more details.

3.1.2.3 Data driven Stochastic Subspace Identification Method (SSI – DATA)

This method makes use of the direct application to the response time series. For more details, reader is referred to [19].

Random Decrement technique is associated with time domain methods but it can also be applied to frequency domain methods. It can be associated with PP, FDD and EFDD to obtain power spectral densities by using FFT algorithm [20]. These methods which are Random Decrement Peak Picking method (RD-PP), Random Decrement Frequency Domain Decomposition (RD-FDD) and Random Decrement Enhanced Frequency Domain Decomposition (RD-EFDD) reduce the noise effect.

A road map to estimate the modal parameters with the help of above mentioned methods

is shown in the figure

Figure 3.1 Block diagram for different methods to estimate modal parameters

Procedure employed in the thesis work is shown in the block diagram in figure 2.3. Road

map to estimate the modal parameters is highlighted by spheres in the figure 3.1. Cross

correlations between response data are obtained by taking the inverse Fourier transform

of the cross PSDs between responses. Cross PSDs are obtained using Welch method

(FFT based method). Modal parameters are then obtained using LSCE method.

4 Parameter identification method

Modal parameter identification method employed in this thesis work is shown as a block diagram in chapter 1. If compared with traditional modal analysis, frequency response functions between input and output are approximated as cross power spectral densities between responses while the impulse response functions are approximated by cross correlations between responses. Prime objective in this chapter is to study the underlying theory related to the estimation of natural frequencies, damping ratios and unscaled mode shapes in the thesis work. Topics under study include power spectra estimation, modal decomposition of PSD, truncation of PSD in OMA, Poles and Residues estimation.

4.1 Power Spectral Density Estimation

Power spectral density is the expected (mean or average) energy density of the signal as a function of frequency.

For zero mean signals, this equals the variance of each frequency sample. In estimating power spectral density (PSD) of a signal, there are two tradeoffs. One is frequency resolution and the other is noise in the signal. During measurements, we have finite length of signal. Two types of errors could occur during PSD estimation. Bias error is concerned with block size of signal while random error is concerned with averages of signal. To obtain a good estimate of PSD, we should have large length of the signal. Only this fact can provide good estimation of PSD. If we take fewer averages in the signal then spectrum will not be smooth. On the other side, if we take small block size, bad frequency resolution could introduce leakage in the spectrum and the necessary information required from the spectrum will be erroneous.

There are two basic methods to estimate the power spectral density explained in sections (4.1.1) and (4.1.2).

4.1.1 Periodogram Method

Simple periodogram method is nothing but quotient of the squared magnitude of the

Fourier transform of the signal and length of the signal. Periodogram method is improved

in Modified Periodogram approach where certain window other than rectangular window

is applied to the signal before taking the Fourier transform. Windows solve the leakage

problem. Bartlett Periodogram averaging introduces the concept of averaging of the

different blocks of the signal. It will decrease variance of the signal at the expense of

resolution. Welch’s Periodogram averaging with overlapping is a method which introduces all the properties of above mentioned methods. This method divides the signal in to blocks and then increases the averaging by taking overlaps of the blocks.

If we have two discrete signals x(n) and y(n) where, (x(n) is considered to be the reference signal) with L samples. We divide it into M blocks with N samples in each block. Consider overlapping of the blocks as well and also let the signals are windowed i.e. multiplied with a window in time domain. The cross power spectra between these two signals is expressed as

∑

=

∗

=

^M

m A

m m w

w

e

yx

C N

X k Y k MB

S

1 ,

( )

,

) 2

( k = 1, 2…….. N/2 +1 (4.1)

m

X

∗w,

is the complex conjugate of the windowed Fourier transform of a block from signal x(n).

C

A

is the window compensation factor given by following relationship,

N n w C

N

n A

∑

⁻

=

= 1 0

) (

(4.2)

B is the equivalent noise band width for the window. It is given by following

e

relationship,

N fs n w

n w N B

N

n N

n

e

.

) (

) (

1 2

0 1 0

2

⎟ ⎠

⎜ ⎞

⎝

= ⎛

∑

∑

−

=

−

=

(4.3)

fs is the sampling frequency of the signal.

w(n) is the discrete window signal N is the block size.

Many windows are available, each one having specific application in signal processing.

Hanning window is used to treat the random signals. So in PSD estimation, we use

hanning window to reduce the leakage effect if the excitation signal is random.

4.1.2 Correlation Method

[23] This method uses the auto or cross correlation between the signals. The correlation function is then Fourier transformed to obtain the corresponding auto or cross spectral densities. Averaging required estimating the PSD is performed in the time domain in this approach. The cross correlation between two discrete signals x(n) and y(n) is given by (x(n) is kept as reference) each with finite block length of N samples,

∑

⁻

−

=

∗

+

= −

^{(1 )}

0

)) ( ) ( 1 (

) (

n N

k

k n x k n y

n N

r (4.4)

With increasing lag n, fewer values are averaged, so more noise is introduced in to the estimated power spectrum. Windowing gives optimum solution to this problem at the expense of less resolution. To cope with it, cross correlation function is multiplied with window and then it is Fourier transformed to find the frequency function. So, final expression for estimation of cross power spectral density is given by following mathematical equation,

( )

∑

−

=

=

^M − M n

jwn

yx

w r n w n e

S ( ) ( ) ( ) (4.5)

4.2 Modal Decomposition of Spectra Densities

As mentioned earlier, we will approximate the cross power spectral densities as the frequency response functions in output only modal analysis. Any structure vibrated with either impulse hammer or shaker can produce a row or a column in FRF matrix. If we measure one row or one column in the FRF matrix, we can synthesize a complete FRF matrix because FRF matrix for linear time invariant system is considered to be symmetric. FRF matrix of size (N m x N m ) is written as

⎥ ⎥

⎥

⎦

⎤

⎢ ⎢

⎢

⎣

⎡

=

) ( )

(

) ( )

( )

(

1

1 11

ω ω

ω ω

ω

j H j

H

j H j

H j

H

m m m

M

N N N

N

L M O M

L

(4.6)

Any ‘N’ degrees of freedom system could be expressed as a sum of ‘N’ single degree of

freedom systems as follows.

[ ] [ ] [ ]

∑

=

×

×

×

+ −

=

^m

−

^{m m m m}

m m

N

r r

N r N

r N r N N

N

j

A j

j A H

1 *

*

)

( ω ω λ ω λ (4.7)

[ H ( j ω ) ] is the FRF matrix, [ ] A

r

is residue and λ is a pole for a particular mode.

r

Poles are complex conjugate pair given by following equation,

2

*

1

,

_{r r r r r}

r

λ ζ ω j ω ζ

λ = − ± − (4.8)

ω is the undamped resonance frequency of mode r, and

r

ζ is the relative damping of

_r

mode r. The residue matrix for mode r, [ ] A can be shown to consist of a product of the

r

mode shape vector { } ψ i.e.

r

[ ] A

r

= Q

r

{ } { } ψ

r

ψ

^Tr

(4.9)

Where Q

r

is a scaling constant that has to be used because the eigenvectors can be arbitrarily scaled (they only tell the shape of the modes).

If equation (4.9) is expanded it gives

[ ]

N r N N

N

N N

r

r

Q

A

⎥ ⎥

⎥ ⎥

⎦

⎤

⎢ ⎢

⎢ ⎢

⎣

⎡

=

ψ ψ ψ

ψ ψ ψ

ψ ψ ψ

ψ ψ ψ

ψ ψ ψ

ψ ψ ψ

2 1

2 2

2 1 2

1 2

1 1

1

L

(4.10)

Equation 4.6 can be written in terms of modal vectors by using equation (4.10) as follows,

(4.11)

The time domain equivalent of the modal model in equation (4.6) is the impulse response matrix [h (t)], where

[ ] ∑ { } { } { } { }

=

×

+ −

=

^Nm

−

r r

T r r r r

T r r r

N

N

j

Q j

j Q H

1 *

*

*

*

)

( ω λ

ψ ψ λ

ω

ψ

ω ψ

[ ] _∑ ( [ ] [ ] )

=

∗

^∗

+

= ^N

^{m r r}

r

t r t

r e A e

A t

h

1

)

( ^{λ λ} (4.12)

As mentioned earlier, input excitation need to be measured to develop FRF. If the excitation forces are considered to be white noise then [7] gives the relationship between the FRF and spectra densities as follows,

[ S

_XX

(ω ) ] = [ H (ω ) ] [ S

_FF

(ω ) ] [ H (ω ) ]

^H

(4.13) Where,

[ S

_FF

(ω ) ] is the spectra density matrix of the (unknown) input forces.

[ S

_XX

(ω ) ] is the spectra density matrix of the responses.

Since the excitation forces are assumed to be white noise sequences, [ S

_FF

(ω ) ] can be considered a constant matrix with respect to the frequency. By substituting equation (4.7) in (4.13) and assuming white noise input, the matrix [ S

_XX

(ω ) ] , with cross power spectra of the responses, evaluated at frequencyω , can be modally decomposed as follows:

[ ] { }{ } { } { } { }{ } { } { }

∑

=

⎟ ⎟

⎠

⎞

⎜ ⎜

⎝

⎛

− + −

− + − + −

=

^Nm

−

r r

H r r r T r r r

H r r r

T r r

XX

j

Q j

Q j

Q j

S Q

1 *

*

*

)

( ω λ

ψ λ

ω ψ λ

ω ψ λ ω

ω ψ (4.14)

Where, { } ψ and

r

{ } Q are respectively the mode shape and operational reference vector

r

for mode r. The reference vector is a complex function of the spectral density matrix of the unknown random input force(s) and the modal parameters of the structure. It should be noted that the modal participation factors Q can not be determined. Since these

_r

factors are required for the normalization of the mode shapes, the operational mode shapes { } ψ remain unscaled i.e. dependent on the unknown excitation forces acting on

r

the structure. Equation (4.14) forms the basis for frequency-domain output-only modal analysis. The order of the model in (4.14) is twice the order of the model in (4.7).

Taking the discrete-time inverse Fourier transform of (4.14) gives the correlation function matrix R(k) for positive and negative time lags k [24]

) (k R =

{ }{ } { }{ } { }{ } { } { }

⎪ ⎪

⎩

⎪ ⎪

⎨

⎧

<

+

≥ +

∑

∑

=

=

0 k for

0 k for

1

* 1

*

*

m

s s r

r m

s r s

r

N

r

H kT r r T k T r r N

r

H kT r r T kT

r r

e Q

e Q

e Q e

Q

λ λ

λ λ

ψ ψ

ψ ψ

(4.15)

It is interesting to note that the casual part of the correlation functions (positive lags) give rise to the first two terms in (4.14) while the non-casual part (negative lags) give rise to the last two terms in (4.14) i.e. (Hermitian conjugate terms). Equation (4.15) is based on the principle of the cross-correlation technique which states that the correlation functions can be expressed as a sum of decaying sinusoids [25]. Each sinusoid has a damped natural frequency and damping ratio that is equal to one of the corresponding structural mode.

4.3 Truncation of the Cross-Spectral Density Matrix

The order of the model in the cross-spectral density matrix (4.14) is twice the order of the model in the FRF matrix (4.7). As a result of this the cross-spectral needs to be clipped, which means that one-half of the length of the cross spectra will be truncated.

For this purpose, cross correlations were computed from the responses. Each cross correlation between responses is truncated i.e. one half of the samples are taken while the remaining part is thrown in accordance to the order of the frequency response function matrix. Then, the truncated Cross correlations are Fourier transformed to obtain power spectral densities whose orders are same as that of frequency response function.

4.4 Poles and Residues Estimation

After acquiring response data, cross correlations and cross PSDs are estimated. Modal parameters are obtained from the poles and residues of the system given in equations (4.8) and (4.9). Curve fitting (Modal parameter estimation) algorithms are used to estimate the modal parameters for each mode of interest in a particular frequency range.

Least squares complex exponential method (LSCE) [17] is used for the estimation of modal parameters in the thesis. This method uses the exponentially decaying impulse responses given by equation (4.12). This method is coded in Matlab toolbox used in the thesis work..

In this method, user is asked to select the frequency range of interest which is used for

curve fitting. Next step in this method is to select the number of poles which are to be

calculated by the estimation algorithm. This assumption is the most decisive in the

algorithm. Usually, the number of poles that give the best result are not exactly the

number of actual poles, but the number of true poles plus a small number of ‘extra’ poles

that take account of errors in the measure FRF’s , sometimes called ‘computational

poles’. In order to select proper number of modes, the algorithm use Stability diagram

and it is obtained using MIF [26]. Examples of stability diagram are shown in the

parameter estimation of a simulation on 4DOF system in figure (6.3) and a bolt channel joint structure in figure (7.5).

After the selection of poles, the driving point FRF is normally used to calculate the residue of that DOF for each mode. This number is then used in the next step, to obtain the complete mode shapes, by applying equation (4.10) [26]. Residues are obtained using least squares frequency domain method in the thesis work.

4.5 Selection of Reference Responses

In output only modal analysis, we are concerned with measuring the responses of the structure. A structure is divided in different points where measurements are to be acquired. Usually, finite element model is used initially so that information regarding nodes is obtained. In traditional modal analysis, either responses are measured with fixed accelerometers (Impulse hammer excitation case) or roving accelerometers (shaker excitation case). In output only modal analysis, measurements are done by keeping one or more locations as reference. This condition depends upon the availability of sensors. If available sensors are enough for the measurement locations then all the experiment can finish in one measurement set. In this case, we can select any location as our reference location and the corresponding cross correlation and cross power spectral density matrices are computed. In case of limited available sensors, references are chosen and other locations are measured in by roving sensors in patches.

Two conditions could occur here. Either the patch of reference responses is moved after two measurements to a new location or they remain fixed at particular location. Former methods are known as jumping or leap frog method while the latter is fixed reference method [10]. A condition applied on leap frog approach is that the reference sensors may land in nodal points of one or more modes of the structure. In the simulation and experimental part of this thesis work, fixed reference method is employed.

Ambient excitations level changes with time i.e. it is of non stationary nature. An extra effort is needed in gluing the different patches of the mode shapes together. Techniques used for time efficient processing of multi-patch output only data under non stationary excitation are not dealt in this thesis work because excitations are considered stationary.

Interested reader is referred to [3] and [27].

5 Normalization of mode shapes

Due to clear advantage of technology and easy testing, output only modal analysis has gained more popularity. As mentioned earlier, to obtain scaled mode shapes we need modal participation factor which is dependent on the excitation. As excitations are not measured so modal participation factors could not be obtained. It concludes that we can get resonance frequencies, damping and unscaled mode shapes. In applications field such as structural health monitoring, structural response simulation and structural modification, scaled mode shapes are demanded so it restricts the usage of output only modal analysis [28].

Some techniques are suggested to get scaled mode shapes. One technique suggested by [29], based on assumption that partition of the inverse of the mass matrix associated with the measured coordinates is diagonal. However, the approach gives answers only when there is a full set of modes, and robustness for a truncated modal space has not been demonstrated. Another approach was suggested by [30] which is based on a more extensive testing procedure that involves repeated testing where mass changes are introduced in the points of the structure where the mode shape is known. This approach seems more feasible as we need information of those modes only which are to be scaled.

Parloo derived an approximate formula for the scaling factor from some basic sensitivity relations in linear dynamics.

Same formula was derived using equation of motion [30]. Theme of the equation is that if we add masses to a structure, it will provide shift in the natural frequencies. So we add masses to the structure at some discrete locations and we get change in resonances for all the modes. Resulting frequency shift for added masses is used in the derivation of formula for the scaling factor.

5.1 DERIVATION

Usually in traditional modal analysis, the modes are scaled such that the modal mass is equal to one, i.e.

{ } ψ

^T

[ ] M { } ψ = 1 (5.1)

Where

{ } ψ is the scaled modal vector [M] is the mass matrix.

It is important to note that the residue matrix [ ] A for a particular pole

r

λ is related to the

_r

modal vectors. For an r-th mode of an N m degree of freedom system:

[ ] A

r

= Q

r

{ } { } ψ

r

ψ

^Tr

(5.2)

[ ]

m r m m

m

m m

r

r

Q

A

⎥ ⎥

⎥ ⎥

⎦

⎤

⎢ ⎢

⎢ ⎢

⎣

⎡

=

ψ ψ ψ

ψ ψ ψ

ψ ψ ψ

ψ ψ ψ

ψ ψ ψ

ψ ψ ψ

2 1

2 2

2 1 2

1 2

1 1

1

L

(5.3)

Where,

Q is the scaling constant that is a function of the scaling of the modal vectors.

r

Using only the q-th column of the residue matrix:

⎥ ⎥

⎥ ⎥

⎥ ⎥

⎥ ⎥

⎦

⎤

⎢ ⎢

⎢ ⎢

⎢ ⎢

⎢ ⎢

⎣

⎡

⋅

= ⋅

⎥ ⎥

⎥ ⎥

⎥ ⎥

⎥ ⎥

⎦

⎤

⎢ ⎢

⎢ ⎢

⎢ ⎢

⎢ ⎢

⎣

⎡

⋅

= ⋅

⎥ ⎥

⎥ ⎥

⎥ ⎥

⎥ ⎥

⎦

⎤

⎢ ⎢

⎢ ⎢

⎢ ⎢

⎢ ⎢

⎣

⎡

⋅

⋅

1 1 1 1 1

1 2

1

ψ ψ ψ ψ

ψ

ψ ψ

ψ ψ

ψ ψ

ψ ψ

qr r

q m

q i

q q

r

mq iq

q q

Q Q

A A A A

(5.4)

For a driving point FRF H

_qq

( ω ) obtained by exciting at point q and measuring the response at point q, the driving point residue for mode r A

_qqr

is given by

2 qr r qr qr r

qqr

Q Q

A = ψ ψ = ψ (5.5)

Equation (4.5) can be used to define modal mass M

r

that is consistent with modal mass defined analytically from the mass matrix.

{ } [ ] { }

r r T

r

M ψ = M

ψ (5.6)

Rather than defining modal mass in terms of an absolute quantity, A

qqr

, modal mass is defined based on Q which is a relative quantity. If the scaling of

_r

{ } ψ

r

is changed, Q will

_r

be altered accordingly. Therefore, the r-th modal mass M of a multi-degree of freedom

_r

system is defined as,

r r

r

j Q

M 2 ω

= 1 (5.7)

r pqr

pr pr

r

j A

M ω

ψ ψ

= 2 (5.8)

Where,

M = Modal mass

r

Q = Modal scaling constant

r

ω = Damped natural frequency

r

From equation (5.7), the modal scaling constant Q

_r

is given by

r r

r

j M

Q 2 ω

= 1 (5.9)

Since, the modes, in traditional modal analysis, are scaled such that the modal mass is equal to one,

r

r

j

Q 2 ω

= 1 (5.10)

Recalling equation (5.5), the driving point residue for mode r, A

_qqr

, is given by

2 qr r qr qr r

qqr

Q Q

A = ψ ψ = ψ (5.11)

Putting the value of Q in equation (5.10) into (5.5), the scaled residues becomes

_r

qr qr

qqr

j

A ψ ψ

ω 2

= 1 (5.12)

The scaled modal vector is { } ψ

r

, is given by

{ } { }

r qr r

r

A

Q ψ

ψ = 1 (5.13)

Usually in output-only modal analysis the mode shapes are just scaled to unity, i.e.

{ } { } φ

^T

φ = 1 (5.14)

{ } { } { }

r ^r

r

A

A

₂

= 1

φ (5.15)

Where,

{ } A

r ₂

is the vector norm of { } A

r

{ } ∑

=

=

^m

i

iqr iqr

r

A A

A

1

*

2

(5.16)

5.1.1 Scaling factor

The scaled and unscaled mode shapes are related by the equation,

{ } { } { } ^φ [ ] { } ^φ

ψ φ

T

M

= (5.17)

Where

{ } ψ is the scaled mode shape,

{ } φ is the unscaled mode shape.

Equation (5.3) implies that the scaling factor is expressed as

{ } φ [ ] { } φ α

T

M

= 1 (5.18)

and

{ } { } ψ = α φ (5.19)

The equation to estimate the scaling factor is derived from the basic equation of motion of a structure subjected to a force { } F ( ) t , i.e.:

[ ] M { } z && + [ ] C { } z & + [ K { } z ] = { F ( ) t } (5.20)

Where,

[ ] M is the mass matrix,

[ ] C is the damping matrix,

[ ] K is the stiffness matrix.

The classical eigenvalue equation in case of no damping is

[ ] { } [ ] { }

1 2

1

1

ω ψ

ψ K

M = (5.21)

Where,

{ } ψ is the mode shape,

1

ω

1

is the natural frequency of any of the modes of the problem related to the mass matrix M.

If a mass change is made so that the mass matrix becomes [ ] [ ] M + Δ M , then the Eigen value equation becomes

[ ] [ ]

( ){ } [ ] { }

2 2

1

2

ω ψ

ψ K

M

M + Δ = (5.22)

Where,

{ } ψ and

2

ω are now the modal parameters of the modified problem.

₂

Subtracting equation (5.8) from (5.7) then gives the equation to approximate in order to obtain an equation for the unknown scaling factor.

[ ] ( { } { } ) [ ] { } [ ] ( { } { }

1 2

)

2 2 2 2

2 2 2 1

1

ω ψ ω ψ ω ψ ψ

ψ − − Δ M = K −

M (5.23)

If it is assumed that the mass change is so small that the mode shape does not change significantly, i.e.

{ } { } { } ψ

₂

≅ ψ

₁

= ψ (5.24)

Hence, equation (4.9) becomes

[ ] { } [ ] { }

2² 2

2 2

1

)

( ω ω ψ ω

ψ M

M − = Δ (5.25)

Pre-multiplying equation (5.25) by { } ψ

^T

results in, { } [ ] { }

2

{ } [ ] { }

2²

2 2

1

)

( ω ω ψ ψ ω

ψ

ψ

^T

M − =

^T

Δ M (5.26)

Considering the orthogonality of the modes,

{ } ψ

^T

[ ] M { } ψ = 1

Hence, equation (5.26) becomes

{ } [ ] { }

2² 2

2 2

1

)

( ω − ω = ψ

^T

Δ M ψ ω (5.27)

Now, recall that in output-only modal testing the modal vector length is scaled to unity

(5.14)

{ } { } φ

^T

φ = 1

Finally, combining equation (4.19) and (4.27) we obtain,

2 2

{ } [ ] { }

2² 2

2

1

)

( ω − ω = α φ

^T

Δ M φ ω (5.28)

And the unknown scaling factor can be obtained as

( )

{ } φ [ ] { } φ ω

ω α ω

T

Δ M

=

₂

−

2

2 2 2

1

(5.29)

When using equation (5.29) to determine a scaling factor, only the mode shape and the natural frequency of that particular mode have to be known. Equation (5.29) gives exact results when the matrix [ ] Δ M is proportional to the mass matrix [ ] M since in this case the modes remain unmodified. [28].

In equation (5.29), both the modified and unmodified mode shapes can be used.

However, the better results are obtained using the unmodified shapes, i.e.

( )

{ }

1

[ ] { }

1 2

2

2 2 2 1

φ φ

ω

ω α ω

T

Δ M

= − (5.30)

or both the unmodified and the modified mode shapes, i.e.:

( )

{ }

1

[ ] { }

2 2

2

2 2 2 1

φ φ

ω

ω α ω

T

Δ M

= − (5.31)

The results of scaling factor obtained with equation (5.30) are independent of the type of

normalization used whereas normalization to the length should be used with equation

(5.31)

5.2 Things to Remember

Mathematical equation for estimating scaling factor shows that we have to add masses at defined locations. We don’t have the actual mass, damping and stiffness values for the defined locations before starting experiment. The added mass should not be very large as it can distort the mode shape. Mass which we add could seem smaller but at the location, at which we add it, may produce different effect as compared to other location. In this way, mode shapes could be erroneous. Adding large masses could cause local damage to the structure. Mass addition also affects stiffness of the structure.

Frequency shifts produced due to mass addition should be significant. Previous research in the field of output only modal testing shows that frequency shifts up to 1% or 2% is considered to be fit to provide good results.

Mass addition should be avoided at node locations. Resulting change in resonance frequency for the particular mode will be very small even for large change in mass. The target is to produce change in resonance frequency with small addition of masses because large change in mass affect mode shapes i.e. mode shapes is deteriorated.

As a conclusion, we have to add masses with small magnitude and should avoid adding them at nodes of the structure.

Multiple mass additions i.e. if we add different masses at more than one location could

improve the results. In this way, mass adding effect distributes among different locations

and we could get reasonable frequency shifts for different modes.

6 Simulation

This chapter deals with application of the theory of Output-Only modal analysis on four degrees of freedom system whose mass, stiffness and damping matrices are known.

System is excited by Gaussian white noise. Responses are simulated and modal parameters are extracted. Moreover, normalization of mode shapes is carried for scaling the mode shapes. Results are compared with the traditional modal analysis.

6.1 Description of System

A four degrees of freedom lightly damped system shown in figure (6.1) is selected for simulation with following parameters.

Table 61 Parameters of the 4DOF system analyzed in simulation

Masses Magnitude

(kg) Damping

Constants

Magnitude

(Ns/m) Spring Stiffness

Magnitude (Nm)

m1 225 c1 600 k1 1500000

m2 175 c2 140 k2 900000

m3 80 c3 175 k3 500000

k4 450000

m4 42 c4 150

k5 200000

Given system is treated as non proportional damped system. Modal parameters extraction for such system with given M, C and K matrices is done by using state space vectors model. As theoretical results used here are just for comparison and the scope of the thesis deal with the methods used in output only modal analysis so the interested reader is referred to Chapter 6 of [32] or Chapter 9 [33].

Theories described in earlier chapters are coded in Matlab tool box. Required results are

interpreted and the figures are plotted for comparing theory and output-only modal

analysis.

Figure 6.1 4 DOF system used in simulation

Given system is excited by Gaussian white noise containing two million samples. Force vector is filtered using fourth order low pass Butterworth filter. Responses are simulated and cross correlations between the responses are evaluated by keeping response at 2 ^nd mass as a reference. In a similar way, cross power spectral densities are computed from cross correlations. Mode indicator function is obtained from cross power spectral densities. MIF is used as a background in the stability diagram for estimation of poles.

Poles are estimated using Complex exponential method. Finally, residues are estimated.

All these methods are coded in matlab tool box.

Modal vectors obtained so far are not scaled. Following mode normalization technique, masses of 4kg and 8kg are added at location DOF1 and DOF4. Same procedure as shown above is followed to obtain the poles and residues for the case of added masses. At this stage, we get poles and residues for with and without addition of masses to the system.

Using these parameters, we get the scaling factor to scale our mode shapes using equation

(5.29).

6.2 Results And Comparisons

Cross correlation between responses at 1 ^st mass and 2 ^nd mass is shown in figure.

Figure 6.2 Cross correlation between Response at DOF 1 and DOF 2

Calculated cross correlations are truncated i.e. half of them are taken because they are of order twice that of impulse response. Stability diagram which shows the location of resonances is shown in figure (6.3).

Figure 6.3 Stability diagram for poles estimation using complex exponential

Estimated poles along with cross power spectra matrices are used to find residues. These residues are unscaled which means that modal vectors are needed to be scaled. Scaling factor is derived by using mode normalization technique defined in chapter 5. This is done by adding small masses of magnitude 4 kg and 8 kg at DOF 1 and DOF 4 respectively. As a result of mass addition, we get a shift in the resonance frequencies is obtained. For the modified structure (structure loaded with masses), we get new poles and residues. Using the evaluated parameters in equation (5.31), we get the scaling factor to scale the residues. Hence, we get scaled FRFs.

Theoretical frequency response functions are calculated using original M, C and K matrices. Comparisons of few modal parameters are shown in figures and tables below.

Figure 6.4 Comparison of Theoretical and Estimated H12

Figure 6.5 Comparison of Theoretical and Estimated H22

Figure 6.6 Comparison of response time histories