FEATURE

by S. O'F. Fahey and J. Pratt

TIME DOMAIN MODAL ESTIMATION TECHNIQUES

I

n a previous technical feature, 2 we described some motives for performing structural dynamic tests and subsequently reviewed some frequency domain techniques for uncovering modal parameters. This time, we take up some applicable time domain algorithms by reviewing the complex exponential algorithm, Pisarenko's harmonic decomposition, Ibrahim's time domain method, and the eigensystem realization algorithm. A numerical example is developed to compare and contrast these techniques.

fers to the unforced response of a system to a prescribed set of initial conditions. It is often denoted simply by the displacement vector y(t). We can model the free decay of a single mode by expressing the homogeneous solution (without forcing) of eq (1) in terms of complex exponentials, or y(t)

=A

exp[~t]

+A*

exp[~*t];

displacement

(2)

where "*" denotes complex conjugate. Substituting eq (2) back into eq (1) and setting ft.t) = 0 yields the characteristic equation

One significant difference between time and frequency domain techniques is the quality of fit as a function of system damping. Time domain techniques generally out-perform frequency domain techniques when the system damping is less than 0.5 percent of critical as a rule of thumb. Frequency domain techniques, on the other hand, generally provide more reasonable results when systems exhibit greater than 4.0 percent damping. This trend can be attributed to how the modal response appears in each domain. For example, a lightly damped response may distribute its effects over the entire duration of a time domain record but produce only a few spectral lines when transformed to the frequency domain.

From the characteristic equation, it is clear that an underdamped (c < 2Vkm) single-mode system, with a nontrivial response (y(t) =1- 0), has two complex conjugate poles

LINEAR MODEL

where

We begin with the time domain description of the structural dynamics of a mechanical system in terms of its free decay and impulse response functions, two common forms of time domain data. Our first assumption is that the structural dynamics in a limited bandwidth can be approximated by a finite number of so-called normal modes. Consistent with this assumption, a dynamic response is taken to be governed by a finite dimensional linear system of second-order ordinary differential equations. Furthermore, this system is described as possessing a principal orientation decoupling individual mode responses such that each mode can be expressed in a single-degree-of-freedom (SDOF) equivalent my(t)

+ cj(t) + ky(t)

= f{t)

(1)

where m, c, k and fare modal mass, damping, stiffness, and excitation; y is a generalized coordinate describing displacement; and the overdot denotes differentiation with respect to time. The impulse response function (IRF) is the system response to a unit impulse. The IRF is mathematically equivalent to the inverse Fourier transform of the frequency response function (FRF) and is commonly denoted h(t). Free decay reS. O'F. Fahey (SEM Member) is an Engineer with Electric Boat Corporation, Groton, CT. J. Pratt (SEM Member) is a National Research Council Research Associate with the National Institute of Standards and Technology, Gaithersburg, MD.

Editor's Note: Time Domain Modal Estimation Techniques is the second in a 2part series on Frequency and Time Domain Methods. If you have any questions or comments about the articles, please contact SEM, [email protected]

[m~ 2

{~" ~n

=

+

+

c~

-2:

±

k}y(t) = 0

(3)

J4~2- ~ (4)

~r =

damping ratio

w, = natural frequency

The first and second time derivatives of eq (2) are sometimes employed to give expressions for free decay in terms of velocity and acceleration respectively. The IRF is the particular solution (with forcing) of eq (1) for the case ft.t) = &(t), where &(t) is the Dirac delta function corresponding to a unit impulsive force at t = 0. The particular solution of eq (1) for this case takes the same form as eq (2), and is called the receptance IRF when it is written in terms of displacement, or the mobility IRF when it is written in terms of velocity. For acceleration, the accelerance IRF takes a different form, namely a Dirac delta function at t = 0 superimposed over complex exponentials. This difficulty is avoided by expressing accelerance in terms of a mobility equivalent. The response of a general multiple-degree-of-freedom (MDOF) system is expressed as the sum of contributing modes, which we write as hpq(t)

=

n

L Apqr exp[~,t] + A;qr exp[~:t]; receptance r=l

(5)

where q, p, and rare reference, response, and mode indices respectively and n is the number of contributing modes. For data originating from an impulse response function, the residual Apqr can be related to modal properties via November/December 1998

EXPERIMENTAL TECHNIQUES

45

TIME DOMAIN MODAL ESTIMATION TECHNIQUES

A

=

pqr

1

2jw, Y1 -

t~

Q pr

qr

(6)

r

where <1>, is the rth mode shape and Q, is called the modal participation factor and is equal to the reciprocal of the modal mass liM,. We point out that the free decay depends not on the forcing but on the initial state of the system. For free decay, one speaks of an initial amplitude Al?'l, which can only loosely be associated with modal participation.

COMPLEX EXPONENTIAL ALGORITHM

= [

~~

~:

YN~2n YN-~n+!

:::

~:::~]

where the subscript serves as our time index, given a constant sample rate t::.t. We have not introduced reference and response indices, as in eq (5), since only one referenceresponse pair is under consideration.

A concern when working with a single reference-response pair is that the dynamics are observable. By observable, we mean that neither response nor reference locations are at a node (zero motion) for any mode of interest, and that all modes of interest have been excited. Having assumed a response dominated by n observable modes, we introduce the characteristic polynomial 2n

2n

k~O

k~!

2: ak(zY' = Il (z -

(8)

zk)

where z denotes a measure of time delay, and the roots zk are equivalent to the system poles

(9) in cycles per second (cps) or Hertz, or (10) in radians per second (for time in seconds), which is our convention unless otherwise noted. We assert that the same characteristic polynomial coefficients relate the columns of the Toeplitz matrix in the case when the order of the model and the system match. Taking a 2n = 1 and partitioning the right hand column, we find the least squares solution for the ak's

... YN-2n+!

Y2n

Y2~+! YN-!

]+{_Y2n+l} { -~~n+2 YN

=

[

exp[X.ftd exp[~ft2 ]

exp[X. 1tN] exp[X.ftN]

:::

exp[X.~t 1 ]] exp[~~t2 ]

... exp[X.~tN]

(13)

reveals our residuals. CEA uses a so-called two-stage regression, requiring the solution of eq (11) then eq (13). Multiple stage-regression is common whenever nonlinear estimation problems can be broken over linearized stages. The problem of linear modal parameter estimation is a nonlinear estimation problem-one reason for having so many applicable techniques. In order to establish consistency of our estimation procedure, the same process is repeated between differing referenceresponse pairs to ensure similar estimates. Systematically, one can employ the least squares complex exponential method of Reference 1 to obtain consistent estimates from single reference multiple response data, or the polyreference time domain algorithm of Reference 9 to obtain consistent estimates from multiple reference multiple response data.

PISARENKO'S HARMONIC DECOMPOSITION Pisarenko's harmonic decomposition (PHD) described in Reference 7 un<;overs an unbiased characteristic polynomial estimate from the covariance matrix (defined below) of a single reference-response pair. The response measure is taken to contain normally distributed noise, e whose mean square value approaches a constant

1

N

2: e~ = N-ooN i~l lim -

1-Lo

so that the measure Yi can be expressed as ui and ei are signal and noise respectively.

(14)

+ E;, where u;

For notational purposes, we express eq (7) as Y = [,y 1y2 ••. Y2n+1l, where the overbar implies a column vector. The covariance matrix is defined, here, as

(15)

~~

CJ.2n-!

EXPERIMENTAL TECHNIQUES November/December 1998

(12)

for which the regressive relation

CJ.o }

(11)

46

exp[X.Iti] exp[~ 1 t 2 ]

=

(7)

Y~

:::

From the pole estimates, we assemble the matrix

A1.N

Our first technique, the complex exponential algorithm (CEA) of Reference 8, is examined through the structure of the Toeplitz matrix for a single reference-response pair

y

The operator (.)+ is called the pseudo-inverse, equal to [(.)T(.W 1(f for real quantities and [(.)B(.)]- 1(.)11 for complex ones, where (.)T, (.)B, and (.)- 1 denote matrix transpose, complex conjugate transpose, and inverse respectively. Having determined these coefficients, we obtain pole estimates using eq (8), and then eqs (9) or (10).

Then, observing

TIME DOMAIN MODAL ESTIMATION TECHNIQUES

E[ .

Y.Yk

l

= {u;u;

+ 110,

uu

' k•

fori = forii:-k

k}

(16)

We pursue a least squares solution between the right and left hand sides of the previous equation, when N - 2 ~ 2n. On the left

where E[.] means expected value, eq (15) becomes

(25)

ufu2n+l -r-

]

u2u2n+l

... urn+lu2:+1

+ 11 (17)

where N > > 1 and 11 approaches (N - 2n)l1o· Examining the structure of eq (17), it can be shown that the eigenvector associated with the minimum eigenvalue of the covariance matrix is an unbiased estimate of the characteristic polynomial coefficients, expressed mathematically as (18)

Finally, poles and residues are calculated according to eqs (9) or (10) and eq (13) as before. Both CEA and PHD uncover 2n observable system poles from the dynamic response of a single reference-response pair. The following sections deal with the premise that each pair of poles implies a system degree of freedom.

IBRAHIM'S TIME DOMAIN Ibrahim's time domain method (lTD) is developed in Reference 3 and is a multiple response algorithm. lTD considers n simultaneous response measures assumed to be dominated by n modes in the form of eq (5). We begin by assembling the response matrices 'fpq(k)

Yp1q, Ym.~k+ll

=

[

Yp.q, ] Yp.q.

Yp2Q,t+t>

YpJq,(k+N-3) Yp2Q,(k+N-3)

(19)

... Yp.q,(k:+N-3)

for which eigenvalue analysis reveals []+[cJ,]{'I'} = [\z\]{'1'}

(26)

and provides frequency, damping and mode shape estimates in a single step! The notation [\. \] implies a diagonal matrix. The mode shapes are the first n entries columnwise in '1', and pole estimates are calculated according to eqs (9) and (10). For a typical modal survey, the number of response DOFs will be greater than the number of system DOFs within a frequency band of interest. The simplest technique to resolve differences between response DOFs and system DOFs is to limit the number of response DOFs considered by analysis. The choice of which response DOFs are considered is arbitrary, but should be representative of the overall modal character. Another technique for resolving differences between response DOFs and system DOFs is principal component analysis as discussed in Reference 5. Principal component analysis is used in multivariate statistics when multiple attributes appear to be related or have similar effects which can not be resolved individually. Here, we speak of a proper orthogonal basis. This orthogonal basis is established using the first n singular vectors of pq(l)T'fpq(l), denoted c from singular value decomposition (a standard matrix computation). The principle components reduce the response DOFs to the system DOFs according to f (k) = 'lpq(k)c. Singular value decomposition is used integr~ly in the following section to resolve between response DOFs and system DOFs.

r

EIGENSYSTEM REALIZATION ALGORITHM

given an excitation, here, associated with a single reference, q and multiple responses, p = {1, 2, ... , n} for k E {1, 2, 3}. Then, two ancillary matrix quantities are formed (transposed from Reference 3) = ['fpq(1)i'fpq(2)]

(20)

= ['fpq(2)j'fpq(3)]

(21)

Each row of either of the proceeding matrices contain two sequential measures of n response measures, which describe the system state (for free decay). Our assumed form of response requires

The eigensystem realization algorithm (ERA) of Reference 4 is a multiple reference multiple response algorithm. We begin our examination of this technique by considering the structure of the response matrices

ypq(k)

Yu. Yt2. ... Ytq.] =

[

Y2~, Y2~,

... Y2q.

Ypt. Yp2.

... Ypq,

...

:

where p, q, and k are response, reference, and time indices respectively; and here ypq denotes an impulse response function. The response matrices can then be COII~pactly expressed as

(22)

= ALN-2"'

(23)

where A is defined by eq (12) and 'It and "' are state representations of the mode shapes. Eliminating A (24)

which is an expression of Ibrahim's method.

(27)

(28) where A, B and C are 2n X 2n, 2n X q and p X 2n matrices. Here, we place the special requirement that A is diagonal. For diagonal A, C and B are called the mode shapes and initial modal amplitudes respectively. The problem of ERA is to construct A, B, and C with the smallest dimension, called a minimum realization, such that eq (28) holds for all k. The size of p and q are the number of response and refNovember/December 1998 EXPERIMENTAL TECHNIQUES

47

TIME DOMAIN MODAL ESTIMATION TECHNIQUES

erence measures. However, the size of 2n need not be related

Finally, mode shapes and initial modal amplitudes are determined by designing two identity matrices

We accomplish our task considering the sum of experimental data in the r X s block matrix (rp X sq matrix)

[Dll2\jlljl-lD-l12pTp] = [ 2n, [QTQD-ll2\jlljl-lDll2] = [ 2n

top nor q.

Hrs(k-

Ypq(k) Ypq(k + 1) [

Ypq(k

+: r-

having used the property that P and Q are self orthogonal = I, QTQ = n, which does not extend to ppT nor QQT. These matrices can be inserted into the right hand side of eq (32) to produce (PTP

1) =

Ypik + s - 1) ] Ypq(k + s)

P[Dll2\jlljl-lD-l12pTp]Dli2A.Dli2[QTQD-ll2\jlljl-lDli2]QT

Ypq(k + r) ::: Ypq(k + s: + r- 2)

1)

from which eqs (32) and (34) and eq (28) and (30) result

also called the generalized Hankel matrix whose entries are matrices themselves, explained by the term block. The choice of the integer values for s and r are arbitrary. However, the smallest dimension of the block Hankel matrix must be equal to or greater than 2n. The authors have found that choosing s equal to the smallest integer greater than or equal to (4n/q) and r equal toN-k - s + 2 is a reasonable, not necessarily optimal, choice. The matrix Hr,(k - 1) can be returned to the response matrix according to

which represent our mode shapes, poles, and initial modal amplitude respectively.

NUMERICAL EXAMPLE We generated a set of numerical data with known properties and added some moderate noise content. This data set involves one reference and five response measures. The dynamics are known to involve three modes, graphically documented by Figure 1 with damping and frequency properties given by Figure 2.

(30) where EJ = [IP, 0P' ... , 0p1pxrp> and IP and 0Pare indentity and null matrices of order p respectively. Similarly, E ~ = [lq,

0q>

••• ,

0q]qxsq·

First, the zero lag block matrix is resolved using singular value decomposition,

1.5

2

2.5

3.5

Time(sj

120

= [P][\D\][Q]T

(31)

where the augmented matrices on the right hand side are rp X sq, sq X sq, and sq X sq from left. to right when rp ~ sq. However, if the dynamics of a limited bandwidth are dominated by n so-called normal modes, then only 2n singular values and vectors are physically significant. By advocacy, the smallest singular values are truncated. One speaks of D, and P and Q as the significant singular values, and left. and right singular vectors respectively. '!'hen, the system eigenvalues can be uncovered by finding A. In according with eqs (28), (30), and (31) (32)

= n-l12pTHrs(l)QD-ll2

ljl-l[D-l12pTHrs(l)QD-ll2]\jl

= [\A\]

EXPERIMENTAL TECHNIQUES November! December 1998

80 40

0

20

40

60

80

100

120

Frequency (Hz)

Fig. 1: Drive Point IRF (top) and FRF (bottom)

AnalFrequency [H~

ping

(33)

[%]

(34)

MAC

and our matrix A is given as [\z 1 ••• z2n \] in the desired diagonal form, where eqs (9) or (10) provide pole estimates in Hertz or radians per second.

48

~

Dam-

where A., is within a rotation, Aljl = ljiA of A; Resolved

A.

100

5'

w 2.80

r•

ytic 34.00 65.30 73.80 1.21 1.09 0.70

CEA

CEA

PHD

lTD

3 34.02 66.28 74.52 1.35 2.79 1.44 0.999 0.921 0.994 0.882

4 34.00 65.40 73.88 1.23 1.38 0.78 1.000 0.999 0.999 0.995

ERA

34.00 65.29 73.76 1.19 1.10 0.71 1.000 0.999 1.000 0.999

34.01 65.31 73.79 1.15 1.07 0.68 1.000 1.000 1.000 Nj_A

34.00 65.31 73.80 1.21 1.09 0.71 1.000 1.000 1.000 0.999

Fig. 2: Parameters and Parameter Estimates

TIME DOMAIN MODAL ESTIMATION TECHNIQUES

Our first three attempts involve estimating frequency and damping properties from the drive point measure shown in Figure 1. The modes of interest are clearly observable in this measure. First, CEA is executed assuming a three mode response. Then, CEA is executed again assuming a four mode response (even though only three mode exist). Finally, PHD is executed assuming a three mode response. Subsequently, residual estimates are obtained for the entire data set using eq (13). Figure 2 demonstrates that the results of our first attempt are poor, mean squared sample correlation (r) is moderate and damping estimates are not close to expectation. In practice, the model order is often over specified due to both uncertainty in model order itself and to compensate for biased model and noise assumptions. So-called computational modes are produced that are not consistent with the physical system. Sometimes pole estimates are plotted versus assumed model order, called a stability diagram, poles estimates that are consistent with increasing model order imply physical significance, whereas unstable poles often imply computation effects. In our second attempt, the computational mode is eliminated and results are improved. Our third attempt using PHD provides good correlation, mode shapes, frequency and damping estimates for this data set. Unbiased estimators provide good results as long as model order can be chosen appropriately. Our last two attempts involve lTD and ERA considering the entire data set. The five response measures are reduced to three measures using principal component analysis and lTD is executed. Frequency, damping and mode shape estimates are obtained. Then, ERA is executed assuming a three mode response. Frequency, damping, mode shapes, and initial amplitude estimates are obtained. Both, lTD and ERA provide good modal estimates for this data set.

CONCLUSION In this brief, we discussed the estimation of modal parameters using time domain algorithms including the complex ex-

ponential algorithm. Pisarenko's harmonic decomposition, Ibrahim's time domain method, and the eigensystem realization algorithm. We emphasize that we have not provided an exhaustive list of methods for uncovering the modal parameters. For instance, References 5 and 6 provide further details of two widely applied methodologies, namely Direct Parameter Identification (DPI) and Auto-Regressive Moving Average (ARMA) respectively. We leave these methods (and others) to the reader with our hope that this brief has served to illuminate key issues facing the experimental modal analyst and provide some avenues to the rich literature on this technologically important subject.

References 1. Brown, D.L., Allemang, R.J., Zimmerman, R., and Mergeay, M., "Parameter Estimation Techniques for Modal Analysis," SAE paper no. 790221 (1979). 2. Fahey, S.O'F., and Pratt, J., "Frequency Domain Modal Estimation Techniques," EXPERIMENTAL TECHNIQUES, 22(5), 3337 (1998).

3. Ibrahim, S.R., and Mikulcik, E.C., "A Method for the Direct Identification of Vibration Parameters from the Free Response," S&V Bulletin, 47(4), 183-198 (1977). 4. Juang, J.-N., and Pappa, R.S., "An Eigensystem Realization Algorithm for Modal Parameter Identification and Model Reduction," AIAA J. Guid., 8(5), 620-627 (1985). 5. Leuridan, J.M., Brown, D.L., and Allemang, R.J., "Time Domain Parameter Identification for Linear Modal Analysis: A Unifying Approach," ASME J. Vib. Acoust., 108, 1-8 (1986). 6. Marple, S.L., Digital Spectral Analysis with Applications, Prentice Hall, Englewood Cliffs, New Jersey (1987). 7. Pisarenko, V.F., "The Retrieval of Harmonics from a Covariance Function,' Geophys. J. R. Astr. Soc., 33, 347-366 (1973). 8. Spitznogle, F.R., and Quazi, A.H., "Representation and Analysis of Time-Limited Signals Using a Complex Exponential Algorithm," J. Acoust. Soc. Am., 47(5), 1150-1155 (1970). 9. Void, H., Kundrat, J., Rocklin, G.T., and Russell, R., ''A MultiInput Modal Estimation Algorithm for Mini- Computers," SAE paper no. 820194 (1982). •

November/December 1998 EXPERIMENTAL TECHNIQUES

49