Identification of modal parameters in vibrating structures

Journal

of Sound and

Vibration

(1988) 125(3), 413-427

IDENTIFICATION

OF MODAL

VIBRATING H.

BARUH

PARAMETERS

IN

STRUCTURES AND

H. P.

KHATRI

Department of Mechanical and Aerospace Engineering, Rutgers University, New Brunsn,ick, New Jersey 08903, U.S.A. (Received 26 January 1987, and in revised form 5 January 1988)

A method is preserited for the identification of the eigensolution associated with vibrating structures. The identification can be performed in the presence of periodic external excitations and noisy measurements. It is assumed that the mass properties of the structure are known, and that the stiffness properties are not known. Linear combinations of the eigensolution obtained with an estimate of the stiffness properties are taken by using unitary transformation matrices to identify the actual eigensolution. The criterion for rotation is based on the number and magnitude of peaks in the frequency spectra of the response of the postulated modal co-ordinates. 1. INTRODUCTION In engineering practice the parameters that describe the evolution of a system are often not accurately known. Sometimes, even the underlying assumptions are incorrect. One has only a postulated model in hand. The task of improving the accuracy of these postulated parameters, based on the system response, is known as parameter identification or system identification. Because of the existence of a large number and different types of dynamical systems, the field of parameter identification has been a very scattered one. In general, an identification method applicable to one type of system is not necessarily applicable to others. Also, there may be several different parameters associated with a system that one may wish to identify. For structures, which are the concern of this paper, it is desirable to identify, among other properties, masses, stiffnesses, boundary conditions, frequencies, eigensolutions, external excitations, damping, non-linearities, etc. The present state of the art in identification theory has been surveyed frequently (see, e.g., references [l-4]). Confining our interest to vibrating systems, we observe that identification methods can be broadly classified into time-domain or frequency-domain categories [4]. Existing work may also be characterized by the property to be identified. In some approaches the mass and stiffness parameters of a structure are identified directly from the system response [5-71. In others, identification of the eigensolution is carried out first [8-111; identification of the mass and stiffness properties is then carried out with use of the identified eigensolution [12-141, or previously known modal information. Identification of external excitations has been considered in reference [15]. This paper is concerned with identifying the eigensolution associated with a vibrating structure. The method introduced here falls into the second category. We propose to identify the eigensolution of a vibrating system by using a combination of time and frequency domain approaches. It is assumed that the mass properties of the structure are known, but that the stiffness properties are not. This assumption is realistic because, while the mass of a structure can be measured with relative ease, the stiffness properties are 413 0022-460X/88/380413+

15 $03.00/O

@ 1988 Academic

Press Limited

414

H. BARUH

AND

H. P. KHATRI

more difficult to assess. The known stiffness properties can be estimated by using the eigensolution identified here in an approach such as the one in reference [12]. Among recently developed methods for identification of eigensolutions are the timedomain methods [8,9], where the free vibration is considered. These methods have been extended to the distributed case [12]. Another approach used for modal parameter identification is that of reference [ll], in which model reduction is also performed. This method, in which only the free response is considered, has been extended to minimize the effects of measurement noise [ 161. In this paper, structures subjected to periodic or non-periodic external excitations are considered. It is assumed that damping is very small. The identification method is applicable to discrete vibrating structures, or to discretized models of spatially continuous systems. This implies that the geometric boundary conditions are also known accurately, since a discretized model is generally constructed by using admissible functions [ 171. It is shown here that the eigenvector matrix of the actual system is related to the modal matrix of the postulated system by a unitary matrix. This suggests that by taking successive unitary transformations of the postulated modes, the actual eigensolution can be identified. The question then arises as to what unitary matrix should be used in each transformation. We use the simplest kind, which is also used in the Jacobi method for computing the eigenvalues of symmetric matrices [ 171, where at each transformation two modal vectors are rotated. To select the rotation angles, the frequency spectra of the modal co-ordinates are considered. If all the system parameters are known accurately the homogeneous response of a modal co-ordinate is sinusoidal, with a frequency spectrum consisting of a single peak. In the event of parameter uncertainties, the modal co-ordinates associated with the erroneous model become linear combinations of the actual ones, leading to frequency spectra with a number of peaks. The rotation angles are chosen so as to minimize the number of peaks in the Fourier transforms of the modes. The modal frequency spectra are generated by using the Fast Fourier Transforms [ 181 of the modal responses, which need to be computed only once, before the identification begins. The identification procedure is applicable to cases where the periodic external excitations are known, such as in a testing environment, or when the amplitudes and frequencies associated with the external inputs are not known. In the latter case, a procedure is outlined to identify the frequencies associated with the periodic excitations. 2. EQUATIONS

OF

MOTION

Consider a discrete system of order n, whose evolution is described by Mf(t)+Kx(t)=F(t)

(1)

where M and K are symmetric and positive definite mass and stiffness matrices of order n x n, respectively, x(t) is the configuration vector of order n, and F(t) is the n-dimensional external excitation vector, which may also include noise. The mass and stiffness matrices can describe a discrete system, or they can represent approximate discretized models of continuous systems such as a lumped parameter or finite element model. Solution of the associated eigenvalue problem consists of real non-negative eigenvalues, related to the natural frequencies by A, = ~3, and corresponding eigenvectors u, (r = 1,2, . . . , n). The eigenvectors are orthogonal and they can be normalized to yield u:Mu, = 6,, ufKu, = 056, (r, s = 1,2,. . . , n). By using the expansion theorem

x(t)=

5 U&(f),

s,(t) =dMx(t),

@a, b)

STRUCTURAL

MODAL

PARAMETER

IDENTIFICATION

415

where q,(t) (r= 1,2,. . . , n) are modal co-ordinates, and the orthogonality the equations of motion can be expressed in the modal form as &(r)+&,(r) where T,(t) = uJF( t) are modal be obtained as the convolution

.., n,

r=l,2,.

= T,(r),

forces. The response sum [17]

of each modal

conditions,

(3) co-ordinate

can then

sin w,( I - T) Tr( T) d7,

(4)

I ql(t)=q,(0)

cos w,t+&(O) sin w,t/w,+-IWYI 0

where ql(0) and 4,(O) are initial conditions. The total response is computed by substituting equation (4) into equation (2a). We will refer to equation (1) as the “equations of motion and to qr( t) as of the actual system”, to u, (I = 1,2, . . . , n) as the “actual eigenvectors”, the “actual modal co-ordinates”. We consider next the case where reliable information is available concerning the mass properties, but the parameters contained in the stiffness matrix are not known accurately. A model for the stiffness matrix is postulated, and this erroneous model is used to analyze the system. With the postulated stiffness matrix denoted by K’, the “postulated equations of motion” become MS(t)+K’x(t)=F(t),

(5)

the physical co-ordinates x(t) being the same in both equations (1) and (5). The associated eigenvectors v, eigensolution is in the form of n eigenvalues A:, and corresponding n), which are orthogonal with respect to M and K’ . We will refer to v, (r=l,2,..., (I = 1,2, . . . , n) as the “postulated eigenvectors”. By using this eigensolution associated with the erroneous parameters the response of the structure can be simulated, by using a relation similar to equation (4), but the simulation results will be incorrect. Our objective is to monitor the response of the system and identify from it the actual eigensolution. The identified eigensolution can then be used to determine unknown parameters in the stiffness matrix or in the actual stiffness distribution [12]. To this end, we introduce “postulated modal co-ordinates”, q:(t) (r = 1,2,. . , n), such that the expansion theorem is expressed as q;(t) = vTMx( t).

x(r) = i v&t), r=l

(6a, b)

Note that these postulated modal co-ordinates q:(f) are not generated from the integration of the equations of motion associated with the erroneous model. They are quantities calculated by using equation (6b), with the actual measurements of the system state. between the actual and postulated eigensolutions. Consider now the relationship Introducing the modal matrices U = [Ii,&

’*utll,

v= [v,vz. **v,],

(7a, b)

one can relate them by

cJ= vc, where C is a square and non-singular transformation unitary. To show this, one can write the orthogonality actual systems with respect to the mass matrix as U’MU

= I,

(8) matrix. It also turns out that C is relations for the postulated and

VTMV = I,

(9a, b)

H.

416

BARUH

AND

H. P. KHATRI

where I is the identity matrix. Substituting equation (8) into equation (9a), and using equation (9b), one obtains C’C = I, so that C-’ = C’ and V = UCT. Similarly, one can relate the actual and postulated modal co-ordinates. Introducing the modal co-ordinate vectors associated with the actual and erroneous models,

* . * %m*9

q(r) = [41(rMr) one can write, using equations

(2) and (6), q’(t) = VTA4x( t).

q(t) = UTA4x( t),

x(t) = Vq’(r),

x(t) = Uq(t),

(10)

q’(t)= [d(tM(t) . * . 4nWlT,

(lla-d) Introducing equation (8) into equation (lld) and using equation that the postulated and actual modal co-ordinates are related by

(9b),

one concludes

(12)

q’(r) = Cq(r).

The aim is to identify the system eigenvalues and eigenvectors (eigenfunctions in the case of continuous systems). Considering the equations above, one concludes that the matrix C needs to be identified. The information available for use in this identification is comprised of the postulated model and its eigensolution, the response vector x(r), and the postulated modal responses q’(t), obtained from equation (lld).

3. THE IDENTIFICATION

PROCEDURE

The identification method proposed in this paper is based on taking linear combinations of the postulated modal vectors (and hence of the postulated modal co-ordinates) in a way so that they successively approach the actual eigenvectors and modal co-ordinates. One first analyzes the response of each mode, which can be expressed as 41(r)=

r=l,2

%(r)+qpr(f),

,...,

n,

(13)

where q,,r(t) = A, cos (wrr - 4,)

(14)

is the homogeneous solution and q,,,.(r) is the particular solution. The total response of each mode is in the form of a single sinusoidal qhr( t), to which excitation terms are added from the particular solution. Note that the excitation may be periodic, in which case the component. The constants A, and 4, response q,(t) will have more than one sinusoidal are functions of the initial conditions and of the excitation. Considering equation (12), one concludes that the postulated modal co-ordinates q:(t) have the form q:(t) = i j-1

c&(r)

= Y? +[A,

cos (Wjr- 4j)’

qpj(t)I,

r=l,2

,*.*,

n,

(15)

I=1

n) denote the entries of C. In other words, the observed time where crj (r,j=1,2 ,..., histories of the postulated modal co-ordinates consist of a summation of sinusoids plus excitation terms. The criterion for identification here will be to transform the postulated modal responses such that the number of sinusoids in q:(t) is reduced to a minimum. We propose to transform the postulated eigensolution into the actual one successively, and at each iteration to use a simple unitary transformation matrix. Because the product of two unitary matrices is also unitary, the aggregate transformation will still be in the form of a unitary matrix. Denoting the initially postulated modal matrix by V = U, , and the matrix C by C = C,, one can write u=

u,c,.

(16)

STRUCTURAL

MODAL

PARAMETER

In each iteration one uses a transformation matrices at the ends of the successive iterations &=

matrix Ri, so that the improved will be

U, = U,R,,

U,R,,

417

IDENTIFICATION

...,

uk+l

=

(17)

UkRk,

where k is the total number of iterations required to identify the actual such that U,,, = U. Considering the first iteration one can write U=

U2Cz= U,RlC2=

modal

eigensolution,

U,C,,

118)

which implies R,C2= The iterative

scheme

Cl,

can be extended

U,,, = U,R, = U,R,R2

. . ’ R,,

CZ = RTC,.

or

(19,20)

to the ith step as

C,,, = R;‘C, = R;R;-,

. . . R;R:C,,

i-1,2

,...,

k,

(21a, b) with U,,, = U and Ck+, = I. The idea behind the approach proposed here is similar to the Jacobi method [17] for computing the eigensolution of a symmetric matrix. The objective here is to diagonalize the matrix C, by successive transformations of the modal co-ordinates. Considering equations (8) and (21b), one concludes that C = C, = R,R2.

. . Rk.

(22)

It turns out that it is not necessary to compute C explicitly, and it suffices to transform successively the modal vectors and modal co-ordinates such that the actual eigensolution is identified. The number of iterations k required to identify the eigensolution is dependent on the accuracy of the postulated model and on the type of unitary transformation at each step. Next, consider the transformation matrices R, (i = 1,2, . . . , k). At each iteration we propose to transform (or to rotate), two modal vectors by an angle. It follows that in order to transform every possible combination of the modal co-ordinate pairs n(n - 1)/2 rotations must be performed. As in the Jacobi method, one may refer to a complete set of rotations of all possible pairs of modal vectors as a sweep [17]. The angle used in a rotation is denoted as I%,.~,where the first subscript describes the number of the iteration and the second and third subscripts indicate the rotated vectors. The entries of R, can then be described as R, = 1,

Rirr = R,,, =

Ri,, = sin 6., ,

COS 6ipyy

j-l,2

,...,

R,,, = -sin

eirr, (23)

n,j#r,j#s,

All other elements of R, are zero. Next the postulated modal vectors can be denoted by v, = ui” and the postulated modal responses by q:(t) = qL’)( t) (r = 1,2,. . , n), where the superscript indicates the iteration number. Considering equations (12), (21a) and (23), one can express the rotated modal vectors and modal co-ordinates at the ith iteration as q!.‘+‘)(t) = -sin

qv+“( t) = cos Oiryq2i’ t)( + sin &,,qY)( t), qy+ ‘I( t) = qJ”( t), Us+”

= COS

6,,u?‘+sin

j=1,2

,...,

n,j#r,j#s,

u(,‘+” = -sin

eiryllj’),

(i+l) _ (8) -u, ” “1

j=1,2

Note that only the rth and sth modal r-s rotation.

vectors

,...,

e,,q’,1’( t) +cos BiT,q:“(f 1,

Olrruv’+ cos 0,,,ut.*),

(24)

n,j#r,j#s.

and modal

co-ordinates

are affected

by an

H.

418

BARUH

AND

H. P. KHATRI

Itremainstodeterminetherotationangle8,,(r=1,2,...,n;s=r+l,r+2,...,n)at each step. Here, we propose to base these rotation angles on the frequency spectra of the postulated modal co-ordinates. To this end, it is useful to examine first the Fourier transforms of the sine and cosine functions. Given a time dependent function h(t), its Fourier transform H(o) is defined by cc H(w)= h(t) emiwrdt, (25) I -uD provided the integral exists for every value of the parameter w. If h( 1) = A cos w,t, then H(w)=(A/2)6(w-w,)+(A/2)6(w+o,), and if h(t)=A sin w,t, then H(w) = i(A/2)6(o - oO) -i(A/2)S(w + w,,). For present purposes it is sufficient to deal with the values of H(w) when w 2 0. Consider next the Fourier transforms (FT) of the actual modal co-ordinates qr(t), which will be denoted by Q(o). Confining attention to the positive values of w, one can write C&(o) = Irn qr( t) e-‘“’ dt=(A, -*

cos &+iA,

sin d,)s(w-w,)+(&,(w),

(26)

where &(w) denotes the Fourier transform of the particular solution, which may have periodic terms. It is assumed that the factor of f is absorbed in A, (r = 1,2,. . . , n). Similarly, the FT of the postulated modal co-ordinates can be expressed as Q?‘(w) = [m q:(t) e-‘“’ dt --a3 =

QQpj(O), j$,Cti(A,COS&+iAj sin~j)s(w-~.;)ci~,

r-l,2

,...,

n,

(27) where it is to be noted that the Fourier transform of the postulated modal co-ordinates may contain many more peaks than Q(o). The criterion proposed for rotation of the modal co-ordinates is based on the knowledge that the FI of the actual modal co-ordinates consists of one peak for the homogeneous solution, and possibly a few more for the particular solution. Assuming that the peaks which correspond to the system frequencies are known or that they can be determined from the frequency spectra of the postulated modal co-ordinates, one can proceed to minimize the number and amplitudes of the peaks associated with the system modes. Upon introducing the terms Re, = GA cos A,

Im, = CA

sin6,

,...,

(28)

r = 1,2, . . . , n.

(29)

r,s=l,2

equation (27) can be expressed as

Q~‘~‘w)=~~~ [(Re~j+ilm~)6(w_wj)+c,jQ,j(w)l,

Note that from the FT of the postulated modes one can measure Rerj and Im,j (r, j = 1,2,..., n) only, without knowing what crj, Aj or 4j are. The FT of the transformed modal co-ordinates becomes Q?+‘)(w) = cos &Ql”(o)+sin Qi+‘ S )(o)

= -sin f3ITS Q (i) * (w)+cos

&Q(si)(~),

0.1TP Q”‘ s (w) 9

i=l,2,..

., k,

(30)

STRUCTURAL

MODAL

PARAMETER

419

IDENTlFICATlON

or QI’+“(W)=d,*6(o-w,)+d,26(o-w*)+...+d,,6(o-w,)+Qb,(w), (31) Q~‘+“(W)=d,~,6(w-w,)+d,*6(w-wz)+...+d,,6(w-w,)+Qb,(w), where

C&(w)

and Q;,<(w) denote

the FT of the corresponding

particular

d,, =

cos

6irsRerr+ sin Bi,,Rec, + i(cos O,,,Im,, + sin @ir,imyr).

d,, =

cos

OirrRers+ sin 8i,xRe,.S+ i(cos O,Jmr~ + sin 9Jm,,

solution

and

), (32)

d,, = -sin

Oi,.SRe,r+

cos

eirsResr+ i( -sin

OJmr,

+

cos Oirslmg,.),

d,, = -sin

OirsRer+

cos

6i,sRe,, + i( -sin

Oirrlmr, +

cos Oir,Jmr,s).

Similar expressions can be developed for other values of d, (i,j = 1,2, . . , n) but for an r-s rotation the interest lies only in d,,, d,,, d,, and d,,. The objective in an r-s rotation is to maximize the amplitude associated with o, in the rth mode, to maximize the amplitude associated with w, in the sth mode, and at the same time to minimize the amplitude associated with w, in the sth mode, and of W, in the rth mode. Mathematically, this can ldsrl and /dJ. Note that be expressed as maximizing ldrrl and Id,,1 and minimizing

l&l # l4.J.

It is easy to show that

so that the sum of the squares of the magnitudes of d,,, d,,, d,, and d,, is independent of the rotation angle &,. This implies that maximizing Jd,,J’+ ld.J2 at the same time minimizes Idrr12+ldyr_12. One therefore chooses the rotation angle to satisfy the following criteria: a(ldrrl’+ Performing

the required

Id.J*)/aeir, algebra

~0,

~*(ldrr12+ 14,l’)/&-s~0.

(34)

gives the value of eirr to satisfy the criteria

tan 2&, = 22/(X

(34) as

- Y),

(35)

where X = Ref,+ Imf,+ Re&+ Imf,, Z = Re,,Re.,, + ImJm,,

Y = Re$ + Im$ + ReSs + 1m:.5, - Re,Re,,

It is clear that in every r-s rotation, the off-diagonal After a sufficient number of iterations d, = 0 (r Z mode becomes free of any contamination from the To demonstrate convergence of the method one P = [P*P2 . . . PJT, P, and P: (r=1,2,..., co-ordinates as

where modal

P, = A,(cos c$, + i sin &),

n) are related

(36)

- ImJm,,

elements, d,, and d,, are minimized. s), so that the response of a certain other modes. can introduce the vectors

P’ = [P’, ZJ; * ’ * lyT,

(37)

to the FT’s of the actual

P: = i c,P, = i (Re, +iim,,), r=, s=,

and postulated

r=1,2

,...,

n, (38)

leading

to the expression P’ = CP.

(39)

H.

420 The norms

BARUH

AND

H. P. KHATRI

of P and P’ are IP12=PHP= $ A;, r=l

so that considering

equations

lP’12=PrHpf=PHCTCp=PHP=

i A;, /:I

(28), (32) and (33) one can write

indicating that the sum of the amplitudes of the peaks in the FT of the system modes is invariant under a unitary transformation. At each iteration C:=l I&]* becomes larger, which makes the corresponding off -diagonal elements & (r, s = 1,2, . . . , n, r f s) smaller, thus leading to convergence. 4. COMPUTATIONAL

CONSIDERATIONS

The Fast Fourier Transform (FFT) is a very useful tool for evaluating the discrete Fourier transform of a function. However, because of the discretization in time and frequency, choice of the sampling period and the number of measurements are very important. One rule of thumb is the Nyquist sampling criterion, which states that the sampling period T should satisfy T < v/w,, [ 181, where wh is the highest frequency of interest. Also, taking a larger number of measurements leads to more accurate results. However, this requires that the system response be measured over a longer time period, which reduces the accuracy of several assumptions made during the mathematical modeling of the structure, such as assuming that the structure is undamped. It is assumed here that the structure is very lightly damped, so that the effects of damping over a short time period are ignorable. Even though the FFT reduces the computational burden associated with evaluating Fourier transforms, the number of required computations is still substantial. Considering equations (30), one can observe that the FFT of the modal coordinates needs to be computed only once, before the identification begins. Once the Q:“(w) (I = 1,2, . . . , n) are calculated, they also are transformed with each rotation, and it is not necessary to compute any more FFT’s. After the identification is completed, one should recompute the frequency spectra of the identified modal co-ordinates to check the accuracy of the identification. There always is some noise associated with a measurement of physical co-ordinates. Instead of x(t) one measures x’(t), where x’(t)=x(t)+w(t)

(42)

in which w(t) denotes measurement noise. In general, w(t) is random and its level of randomness is related to the magnitude of x(t). The question remains as to the nature of the frequency spectrum of w(t). Because w(t) is random its plots are very wrinkled functions, resembling a summation of a large number of sinusoids, all with a relatively uniform amplitude. It follows that each peak in the Fourier transform of the modal co-ordinates is affected by about the same amount. One intuitively expects that a moderate amount of noise should not affect the identification results. The identification procedure was described in equations (27-36) in terms of displacement measurements for every degree of freedom. Examination of these equations indicates that the identification can be carried out by velocity and acceleration measurements as well. Also, if M and K are derived from very high order finite element models, it becomes unfeasible to rotate all possible modal vector pairs. In such cases, only the lower modes should be transformed.

STRUCTURAL

MODAL

PARAMETER

IDENTIFICATION

421

In many cases the amplitudes and frequencies of the external excitations acting on a structure may not be known. If the excitation is periodic, the frequency spectra of the modal co-ordinates will contain peaks associated with the excitation, which need to be separated from the peaks belonging to the modal quantities. To achieve this, consider the postulated model. If the postulated model is realistic, the nature of its eigensolution can be assessed. Mathematically, one can describe this as the matrix C having larger terms in its diagonal elements than in its off-diagonals. One would then expect the first peak to be the largest one in the frequency spectrum of the first postulated mode, the second peak to be the largest in the frequency spectrum of the second postulated modal co-ordinate, and so on. Use of this knowledge helps identify the frequencies associated with the excitation, as will be illustrated in the next section. Note that if a wrong peak is identified as belonging to the external excitation, the identification procedure will not converge. The rotation angle may become infinitesimally small, but the number of peaks in the frequency spectra of the modes will not be minimized. This is the main difference between the approach proposed here and the Jacobi method. A similar statement can be made for the case when some of the underlying assumptions associated with the mathematical model are incorrect, such as when the mass parameters are not known accurately. These problems associated with the implementation of the procedure have been discussed in detail in reference [19]. Considering the rotation of the postulated modal co-ordinates, one can note that at each sweep, ( n2 - n)/2 transformations are performed, which permits rotations between every possible modal co-ordinate pair. The iterations are continued until all angles 0,,, are smaller than a threshold value in a certain sweep or until a maximum number of sweeps is exceeded. The reason a limit on the number of sweeps is chosen is to avoid possible errors associated with recording the peaks in the Fourier transforms, which may prevent convergence. It should also be noted that criteria other than the one described here for rotation of the postulated modes can be used to implement the identification procedure. The identification procedure can be summarized as follows: (1) develop a postulated model, solve the associated eigenvalue problem, and obtain the postulated modal vectors n; (2) record the displacement (or velocity or acceleration) measurements; v,,r=1,2 ,..., these quantities can be recorded directly if an actual experiment is conducted, or they can be obtained by simulation of the actual structure on a digital computer; (3) compute the postulated modal co-ordinates qL( t) by q:(t) = v~Mx( t) (r = 1,2, . . . , n); (4) obtain the Fourier transforms Q:(w) of q:(t) (r = 1,2, . . . , n) and record the peaks; the Fourier transforms are computed only once, before the identification begins; (5) identify and separate in Q:(w) peaks associated with the external excitations; (6) rotate the postulated modes until all rotation angles &, are less than a threshold value or until a maximum number of sweeps is exceeded.

5. ILLUSTRATIVE

EXAMPLES

Two models of vibrating systems can be used to illustrate the identification procedure. One is a three degree-of-freedom model and the other is an eight d.o.f. mass-spring system. The reason a low-order example is also presented is because the actual and identified modal vectors can be tabulated and compared with more ease. The models are shown in Figures 1 and 2. For the first case, the following parameters were chosen: mass coefficients, M, = 1 .O, K, = 3.5, K, = 6.0, K3 = 4.0; postulated MZ = 2.0, M3 = 1.5; actual spring constants, springs, K, = 6.5, K2 = 3.0, K3 = 4.0. A forcing of F(t) = 1.0 sin t was added to the first

422

H. BARUH

AND

H.

P. KHATRI

Figure 1. The eight d.o.f. system. Actual system: M, = 5, Mz = 3, M3 = 4, M4 = 5, M, = 3, M,, = 4, M, = 5, M, = 3, K, = 100, K, = 80, K, = 50,K, = 200, K, = 30, K, = 50;postulated data: K i= 80, Ki = 60, KG = 100, K;=200, K;=50, Kk=20.

Figure 2. The three d.o.f. system. Actual system: M, = 1.0, M, = 2.0, Mx = 1.5, K, = 3.5, K, = 6.0, K, = 4.0; postulated data: K i= 6.5, Ki = 3.0, K; = 4.0.

mass for the forced input cases. For the eight d.o.f. model, the system parameters were chosen as follows: mass coefficients, M, = 5, M2= 3, M,=4, M4= 5, M,=3, M6=4, M, = 5, M8 = 3; actual spring constants, K1 = 100, Kz = 80, K, = 50, K4 = 200, K, = 30, K6 = 50; postulated springs, K, = 80, Kz = 60, K3 = 100, K4 = 200, K5 = 50, K6 = 20. For the forced vibration cases a force F,(t) = sin t was applied to mass 1 and a force F2( t) = 0.5 sin 2t was applied to mass 2. In the free vibration case, the response obtained was that for given initial displacements. As can be seen, the initial estimates of the spring constants are very inaccurate. To simulate the effects of measurement error, noise was added to the displacement measurements in the form x;(t)=xj(t)+AjRj(t),

j = 1,2, . . . ) n,

(43)

where Xj( t) is the actual displacement of the jth mass, xJ( t) is the measured displacement, Aj is a noise factor, and Rj is a random variable with a uniform distribution of [-0*5,0*5]. The parameter A, was chosen as the product of the maximum displacement amplitude of the jth sensor and a selected noise level (expressed as a percentage). For example, if during a certain cycle the maximum displacement amplitude is 0.3, and 10% noise is considered, A, is taken as 0.3 x 0.1. Even though measurement noise is generally treated as Gaussian, it was felt here that the model assumed above is more realistic, because the magnitude of the noise is related to the amplitude of the displacement. During system simulations a sampling period of O-122 s was used for the three d.o.f. case, and a period of 0.024 s was used for the eight d.o.f. case to generate the FFT plots. The sampling periods were chosen so as to satisfy the Nyquist sampling rate criterion [18]. The threshold value for the rotations was taken as 10m6 rad, and the maximum number of sweeps was chosen as six. The natural frequencies of both models were identified with good accuracy from the frequency spectra, and they are given in Tables 1 and 2. In order to evaluate the accuracy

STRUCTURAL

MODAL

PARAMETER TABLE

IDENTIFICATION

423

1

Exact and identijied frequencies for three d.o._tYcase Mode no.

Exact

Identified

1

0.71850

2

2.10784

2.11115

3

3.49391

3.49345

TABLE

0.71628

2

Exact and identljied frequencies for eight d.0.f: case Mode no.

Exact

Identified

2.21867 3.79990 5.42702 6.61996 7.63204 9.45154 11.79351 14.42000

2.23053 3.80132 5.43496 6.59734 7.63407 9.45619 11.78097 14.41991

of the identified eigensolution, one can define a modal error coefficient e, as e, = [Ill, -u;k+r)ll”2, where ej is a measure of the error in identifying the jth mode and where uj is the exactjth eigenvector and u:~+~) is the identifiedjth eigenvector. Table 3 compares the exact, postulated, and identified eigenvectors, as well as the identification error for the three d.o.f. model. Various levels of measurement noise and excitation were considered. One can observe that moderate levels of noise and presence of a periodic external excitation have negligible effect on the accuracy of the identification. Actually, identification results seem to be better for the case of periodic excitations. In some cases there seems to be a very small degradation in the accuracy as the noise level is increased, which would be expected. However, in other cases the accuracy of identification did not change as the noise level was increased. The rates at which the eigensolutions were identified were about the same for different noise levels and excitations. Naturally, a good initial estimate speeds up convergence. No major differences in convergence were observed between the periodic excitation and the free response cases. The slight differences between the identified and actual eigenvectors are due to errors in recording peaks in the FFT plots. This can be corrected by recomputing the FFT’s of the identified modes, treating the identified modes as comprising the new postulated model, and by conducting the identification again. However, inaccuracies involved in the identification were too small to necessitate such an analysis. Next, for the eight d.o.f. case, modal errors and the actual and identified eigenvectors for different cases of noise and excitation can be compared. The modal error coefficients are given in Table 4. As can be seen, the results are almost identical to the results obtained for the three d.o.f. case, corroborating earlier conclusions regarding the accuracy of identification, effects of noise, and presence of external excitations. Here, however, the effects of increasing the noise level are more visible.

424

H.

AND

BARUH

H. P. KHATRI

TABLE

3

Comparison of actual, postulated and identiJied eigenvectars and error norm for the three d.o.J: case Actual eigenvectors 0.3063536 0.4587008 0.5688202 Original

postulated

0*1690614 0.5045650 0.5551255

0.4553880 0.3838158 -0.5761902

0.8359242 -0.3771989 0.1054282

eigenvectors 0.8342970 0.1864381 -0.3953312

O-5247539 -0.4589717 04496838

Case 1: no forcing, 1% noise Identified eigenvectors 0.3062543 0.4584383 0.5691379 e, = 4.2381 x 10p4;

0.4556438 0.3840114 -0-5758815 e2 = 4.4603 x 10e4;

0.8358212 -0.3773189 0~1054004 e3 = 1.6049 x 1O-4

Case 2: no forcing, 5% noise Identified eigenvectors 0.3061570 0.4584846’ 0.5691231 e, = 4.2087 x 10e4;

0.4555517 0.3840506 -0.5758953 e2 = 4.1099 x 10m4;

Case 3: no forcing, 10% Identified eigenvectors 0.3060350 0.4585425 o-5691047 e, = 4.5547 x 10e4;

noise

0.4554354 0.3841002 -0.5759124 e2 = 4.0030 x 10w4;

Case 4: Sinusoidal forcing, Identified eigenvectors 0.3065806 0.4585898 0.5688580 e, = 2.5554~

10m4;

0.3063935 0.4585502 0.5689677 e, = 2.1463 x 10-4;

5%

10%

0.8358094 -0.3773517 0.2053061 e3 = 2.2676 x 10e4

noise

O-4553265 0.3839003 -0.5761475 ez= 1.1297~ 10d4;

Case 6: Sinusoidal forcing, Identified eigenvectors

0.8360151 -0.3771019 0.1054108 e3 = 1.3391 x low4

1% noise

0.4554459 0.3837983 -0.5761752 e2 = 6.2256 x 10e5;

Case 5: Sinusoidal forcing, Identified eigenvectors 0.3065089 0.4585746 0.5689001 e,=2.1553~10-~;

0.8359070 -0.3772228 0.1054053 e3 = 3.7067 x lo-’

0.8359007 -0.3772663 0.1052306 e,=2.1022x10-4

noise

0.4551324 0.3840653 -0.5761031 e, = 3.6746 x 10p4;

0.8360487 -0.3771280 0.1051077 e3 = 3.5489 x 10m4

STRUCTURAL

MODAL

PARAMETER

425

IDENTIFI(‘ATION

4

TABLE

Comparison of error norms for eight d.o.j case No forcing, no noise

Mode number

1 2

1.3151 7.7269 1.3688 1.1751 1.1548 5.8386 6.0053 3.9356

3 4 5 6 7 8

x x x x x x x x

No forcing,

Forcing, no noise

5% noise

1K2 10-j 10m3 lo-? lo-’ 1O-4 10m4 10m4

1.5628 9.7615 9.0713 1.2756 I.2905 8.4800 1.2221 6.7768

x x x x x x x x

lo-’ 1O-4 10-j 1O-3 10m3 1O-4 lo-’ 10m4

1.2483 x 7.4302 x 1.3767 x 1-1794x 1.1553 x 5.6429 x 6.0554x 3.8808 x

Forcing, 5% noise

10-j lO--4 10m3 1om3 1om3 10 -4 1O-4 lo-’

1.4717 x 9.8849 x 1.3777 x 1.4152x 1.5165 x 9.0165 x 8.0700 x 5.7344x

IO i 10 ’ to-3 10~’ IO-” lo-~” 10 -’ 10 ‘I

can consider identification of peaks associated with external excitations whose frequencies are not known. In Figures 3-5 plots of the frequency spectra of the postulated modal co-ordinates for the three d.o.f. case are shown. The frequency of the sinusoidal excitation is between the first and second natural frequencies. Considering the three plots, one can observe that the amplitude associated with the jth natural frequency is highest in the frequency spectrum of the jth mode (j = 1,2,3), which is to be expected. One can then examine each figure individually, assuming that one does not know the Finally,

one

0.5

0

I

I.5

2

2.5

3

3.5

(Rad/s)

Figure

0.125

3. Frequency

spectrum

of first postulated

mode

-

0.1 2 .g

0.075

-

0.05

-

0.025

Or‘ 0

J\ 0.5

I

1.5

2 w

Figure

4. Frequency

*

&

6

spectrum

2.5

3

3.5

(Rad/s)

of second

postulated

mode.

4

H. BARUH

426

AND

H. P. KHATRI

0.06 s .10 0 0.04

I

I.5

2

2.5

3

3.5

4

(Rod/s)

Figure

5. Frequency

spectrum

of third postulated

mode.

frequency of the external excitation. Considering Figure 3 alone does not help determine which peak belongs to the sinusoidal input. Figure 4 indicates that because the third peak has the highest amplitude, one of the first two peaks must be associated with an external excitation. Figure 5 leads to the same conclusion. By reconsidering Figure 3, one can then realize that the first peak belongs to mode one, because the amplitude of the second peak varies erratically for different mode numbers, which is observed from Figures 3-5. Amplitudes associated with periodic excitations decrease uniformly as the mode number increases for cases of w > w,, where w is the driving frequency. In any event, as stated before, if a wrong peak is identified as belonging to the excitation the identification procedure does not converge. 6. CONCLUSIONS A method has been presented for the identification of eigensolutions of vibrating systems. It is assumed that the mass properties are known, and the stiffness properties are not known. The identification is based on unitary transformations of the postulated modal co-ordinates. The criterion for rotation is based on the number and amplitude of peaks in the frequency spectra of the postulated modal co-ordinates. The identification can be implemented in the presence of known or unknown periodic excitations and noisy measurements, and it is found to be relatively insensitive to substantial amounts of noise. 7. ACKNOWLEDGMENT This work was supported in part by the U.S. National Science Foundation, MEA-85-01877.

Grant No.

REFERENCES 1. A. BERMAN

1979 Shock and Vibration Digest 11(l). Parameter identification techniques for vibrating structures. 2. V. STREJC 1981 Automatica 17(l), 7-21. Trends in identification. 3. C. S. KUBRUSLY 1977 International Journal of Control 26(4), 509-535. Distributed parameter identification, a survey. 4. L. LJUNG and K. GLOVER 1981 Automatica 17,71-86. Frequency domain versus time domain methods in system identification. 5. H. T. BANKS and J. M. CROWLEY 1983 Proceedings ofthe 1983 American Control Conference, June 22-24, 1983, San Francisco. Parameter identification in continuum models. 6. L. A. BERGMAN, A. L. HALE and J. C. GOODING 1985 American Institute of Aeronautics and Astronautics Journal 23, 1234-1235. Identification of linear systems by Poisson moment functionals in the presence of noise.

STRUCTURAL

MODAL

PARAMETER IDENTIFICATION

427

I. D. C. SAHA and G. PRASADA RAO 1980 Proceedings of the Institution of Electrical Engineers 127, 44-50. Identification of distributed systems via multidimensional distributions. 8. S. R. IBRAHIM 1983 American Institute of Aeronautics and Astronautics Journal 21, 446-451. Computation of normal modes from identified complex modes. 9. S. R. IBRAHIM and E. C. MIKULCIK 1977 Shock and Vibration Bulletin 47, 183-198. The experimental determination of vibration parameters from time responses. 10. D. J. EWINS and P. T. GLEESON 1982 Journal of Sound and Vibration 84, 57-79. A method for modal identification in lightly damped structures. 11. J-N. JUANG and R. S. PAPPA 1985 Journal of Guidance, Control, and Dynamics 9, 620-621. An eigensystem realization algorithm for modal parameter identification. 12. H. BARUH and L. MEIROVITCH 1985 Journal ofSound and Vibration 101, 551-564. Parameter identification in distributed systems. 13. M. BARUCH 1982 American Institute of Aeronautics and Astronautics Journal 20. 441-442. Correction of stiffness matrix using vibration tests. 14. A. BERMAN 1979 American Institute of Aeronautics and Astronautics Journal 17, 1147-l 148. Mass matrix correction using an incomplete set of measured modes. 15. H. BARUH and L. SILVERBERG 1986 JournalofSoundand Vibration 108,247-260. Identification of external excitations in self-adjoint systems using modal filters. 16. J-N. JUANG and R. S. PAPPA 1986 Journal of Guidance, Control and Dynamics 9, 294-303. Effect of noise on modal parameters identified by the eigensystem realization algorithm. 17. L. MEIROVITCH 1980 Computational Methods in Structural Dynamics. Sijthoff and Noordhotf. 18. E. 0. BRIGHAM 1974 The Fast Fourier Transform. Englewood Cliffs, N.J. Prentice Hall, Inc. 19. H. BARUH and J. BOKA 1988 Proceedings of the 1988 AIAA Structural Dynamics Conference, April 18-20, 1988, Williamsburg, VA. 1442-1456. Implementation problems associated with modal identification in flexible structures.