Analytic evaluation of spectral moments

Journal of Sound and Vibration (1988) 124(3),479-488

ANALYTIC EVALUATION OF SPECTRAL MOMENTS M.

DI

PAOLA AND

G.

MUSCOLINO

Dipartirnento di Ingegneria Strutturale e Oeotecnica, Unioersita di Palermo, Viale delle Scienze, 90128 Palermo, Italy (Received 30 October 1986, and in reoised form 17 November 1987) In this paper an analytic procedure that drastically reduces the computational effort in evaluating the spectral moments of the response of multi-degree-of-freedom systems is presented. It is shown that the cross-spectral moments of any order of two oscillators subjected to a filtered stochastic process can be obtained in a recursive manner as a linear combination of the spectral moment of each oscillator up to the third order separately taken. A numerical procedure is also presented in order to evaluate such first few spectral moments. 1. INTRODUCTION

Many forms of dynamic excitation, such as waves, wind, etc., can be characterized as stochastic processes; it follows that the response of the structures loaded by such agencies can be described in a probabilistic sense. In this framework, structural safety can be evaluated by using the spectral moments (SM) [1]. The latter represent the moments of the one-sided power spectral density (PSD) function with respect to the frequency origin. For multi-degree-of-freedom (MDOF) linear systems, by using modal analysis, the SM of nodal variables can be obtained as a linear combination of the modal SM. The latter can be computed by integrating, in the frequency range, the cross-P8D (CPSD) of the modal output, given as the product of the excitation PSD by the transfer functions of two structural modes. The modal 8M are given in a closed form solution only in some particular cases [2,3], while usually they are evaluated numerically. Moreover, to obtain the nodal statistics of the output in general, m(r+ 1)(m + 1)/2 cross-modal SM must be computed (r being the order of SM required and m the lower natural modes selected for the analysis). It follows that a very large computer time is necessary, due to both the high number of modal 8M required and to the dense frequency discretization in the frequency domain. The aim of this paper is drastically to reduce the above-mentioned drawbacks, which often restrict the use of random vibration theory. Starting from the definition of the one-sided CPSD and applying the Wiener-Kinchine theorem twice, the CPSD is again obtained in a different but equivalent expression, namely as a linear combination of the squared modulus of each modal oscillator in decoupled form. Furthermore, the decoupled form remains unchanged when one takes the Fourier transform of the time derivatives of the cross-correlation function, and the coefficients of the squared modulus of each modal oscillator are obtained in a recursive manner [4]. It follows that the cross-modal SM of any order are here obtained in a recursive form as a linear contribution of the SM up to the third of each modal oscillator order separately taken. In this way, to obtain the nodal SM, only 4m modal 8M need be evaluated by numerical integration. It will be emphasized that the aforementioned properties are independent of the shape of the PSD of the input. 479 0022-460X/ 88/150479+ 10 $03.00/0

© 19&8 Academic Press Limited

480

M. Dl PAOLA AND G. MUSCOLINO

At the very least, the dense frequency discretization for the evaluation of the first four 8M having equal indices can be avoided by using the Pfaffinger method [5], improved by the use of the recursive formulas presented here. By using the procedure outlined above, a large number of error sources are eliminated. The cross-8M of any order can be obtained exactly, simply by knowing the 8M of the single oscillators up to the third order, and these quantities can be evaluated in a very accurate manner by using third order polynomials to describe the PSD of the input.

2. PRELIMINARY CONCEPTS The equation of motion of an n-degree of freedom linear structural system subjected to a forcing function f( t) is written in the form:

Mx+ ex + Kx = LJ(t),

(1)

where M, C and K are mass, damping and stiffness matrices respectively, x is the displacement vector, the upper dot means time differentiation, L is an influence vector and J(t) is a zero mean stochastic process. For classically damped systems the usual modal transformation (2)

q = cIlx, where cIl is the modal matrix normalized with respect to M, yields

q+ 4.q+n 2 q = vf(t).

(3 )

Here, A and n are diagonal matrices listing 2l:j W i and Wj respectively, l:i being the percentage of critical damping and Wj the natural radian frequency. The vector v is given as

(4) The ith equation of the system (3) can be written as

iij+2l:jwjQj+WTQj = vJ(t).

(5)

If f( t) is a stationary stochastic process the modal vector q (t) and the nodal vector x( t) are stationary stochastic processes too, and both are characterized by the power spectral density function (PSD). The PSD of the input and of the output are defined as

E[F*(w,T)F(w,t)] S ()-l' E[Qt(w,T)Qs(w,T)] = 1' S( 'Jf ta ) im 7T T 'qkq, W - T_c(I im 7T T_oc T

,

(6,7)

where the asterisk denotes the complex conjugate, and F(w, T), and Qj(w, T) (i = k, s) are the truncated Fourier transform of the processes f(t) and qj(t) respectively, i.e.,

F(w, T) =

LTT f(t) e-

iw

,

dt,

Qj(w, T) =

LTT q;(t) e-

iw

,

dt,

(8)

i being the imaginary unit. It is worth noting that the PSD defined in equations (6) and (7) are defined in the range -00 < w < 00, It is widely understood that in order to evaluate some quantities of engineering interest, such as first excursion probability, threshold crossing, etc., the so-called one-sided PSD must be introduced. The latter, which from here on will be denoted by the symbol G(w), is related to the corresponding PSD defined in equations (6) and (7) through the following relationships:

Gff{w) = 2U(w)Sff(W),

Gqkq,(w)=2U(w)Sqkq,CW).

(9,10)

481

EVALUATION OF SM

U( w) is the unit step function. The moments of the one-sided PSD about the frequency origin are the spectral moments (8M): i.e., A."ff== foe, w'2U(w)Sff(W) dw -00

A."qkq,(W) =

foo

= rOD w'Off(W) dw,

(11)

Jo

w'2U(w)S'M,(W) dw =

foo W'Oqkq,(W) dw.

(12)

0

-00

Due to the fact that the SM of the vector of the nodal variable x can be computed as linear combinations of the SM of modal variables q in the form m

,n

L L cPjk cPjSA. r,qkq"

A.r.xjxj ==

(13)

k= 1,=1

attention will hereafter be focused on the evaluation of the SM of modal variables. The CPSD defined in equation (10) is related to the PSD of the input through the relationship 0M,(W) = vkv.ijqkq,(W)Sff(W),

(14)

where Vj (j == k, s) is the jth entry of the vector v, and Oqkq,(W) = 2U(w )Ht(w )H..(w),

(15)

Hj(w) (j = k, s), being the transfer function: i.e., H/w)== 1/[(wJ-w2)+2i~JwJw].

(16)

The function Oqkq,(W) defined in equation (15) can be regarded as the CPSD of the two oscillators qk(t) and qs(t) subjected to a white noise (WN) input having a one-sided PSD. Inserting equation (14) in equation (12), one obtains A.r,qkq,=VkV,

{'C' W'Oqkq,(w)Sff(w)dw,

(17)

and due to the fact that in equation (15) the product between Ht(w) and Hs(w) appears in a coupled form, a large number of numerical integrals (of order (r+ l)m(m + 1)/2, m being the number of lower modes selected for the analysis), must be evaluated, with dense discretization in order to obtain good accuracy. However, in the next section the compact form given in equation (15) will be written in a decoupled form involving the transfer functions of each oscillator separately taken. 3. DECOUPLING OF CROSS-PSD

The cross-correlation function corresponding to the one-sided PSD Gqkq,(w) defined in equation (15) can be obtained, according to the Wiener-Kinchine theorem, in the form ....(0)

Rqkq,(T) =

"'" Joroo Oqkq,(W) e

iWT

dw,

(18)

where the superscript in parenthesis means derivative with respect to the time T. The correlation function defined in equation (18) is a complex one; the real part has the form Re [R~~~, (T)] = {'IT exp (-~sWsT)[ ~;'kS cos.WsT + ~~,kS ~in _WsT], , 'lTexp«(kWkT)[ao,skcoswkT-{3o,sksmwkT],

T;;;' OJ, T=;;;O

(19)

where Wk and Ws are the damped radian frequencies Wk = w kJl - ~~,

Ws = wsJl - ~~.

(20)

482

M. 01 PAOLA AND G. MUSCOLINO

The other quantities are as follows:

et; = 2(wi -

a;'ks = 4«(kWk + ("w.,)/ Kks>

w; +2?k'"WkW" + 2(;w.~)/(w.,Kk..)'

K k., = (wi - w;)2+ 4(k?s(wk + W;)WkWs + 4('~ + ?;)w~w;.

(21)

aO:sk, f3~sk are evaluated by changing the indices in equation (21) (the superscripts plus and minus may be omitted in equation (21) for k ¥ s, but for k = s the superscripts are useful to avoid confusion). The imaginary part of R~~~,Cr), defined in equation (I8), is the Hilbert transform of the real part defined in equation (19) [6,7]: i.e.,

f'"

I Im [R(O) (r)] = HT {Re [R(O) (r)]} = q.q,

q.q,

7T'

c-eo

Re [R(O) ( )] M, Pdp. T -

p

(22)

The symbol HT { . } denotes the Hilbert transform operator. Furthermore, due to the presence of the factor WrOqkq,(w) in the integral given in equation (17), it is convenient to form the derivatives of order r of the correlation function defined in equation (I8). Equation (19) shows that the cross-correlation function is defined through different functions in T> 0 and r < O. It follows that the time derivatives of order r must be written in the form 7T exp

r~4

r} ( )] _ R e [R- (q.q., r -

(-("w,r)[ a:

t: cos W"T + f3;'k.r sin w"r],

(

d

{

w,-m-I,k"

m8(r»)

_

+

m , 0 <1'<0 d1' . 7T exp «(kwk1') [a ~sk cos WkT - f3 ~sk sin Wkr], where B(r) is the Dirac delta function and

7T I...

(23)

m=O

f3~k."i= -Cs w sf3 :- 1,ks- WsO:; -

l .ks ,

a;'ks= -?wSa:-l,kS+ ws/3 ;- I.ks> +

-

(24)

wr-m-I.ks = a r-m-l,ks - a ,-m-I,sk'

Notice that a;'k,,=-a~.'k and f3;'ks=-f3~sk for n odd and k=s, while a:'kS=a~Sk and f3;'ks = f3 ;'sk for n even and k :=: s. Equation (23) makes it possible to compute all derivatives of the correlation function in a recursive manner. In order to explain the presence of the Dirac delta functions in equation (23), in Figure I the cross-correlation function and its time derivatives up to the fourth order are depicted. These figures show the behaviour of the cross-correlation function of the given system subjected to a real white noise input. The cross-correlation function and its first derivative (Figures I(a) and I(b» are continuous functions at T = 0, the second derivative (Figure l(c» exhibits a discontinuity of the first kind at T = 0, while the thi rd derivative (Figure 1(d) exhibits a jump and therefore the fourth derivative (Figure l(c» exhibits a jump and a Dirac delta function at r= O. It follows that, with increasing order of cross-correlation differentiation, singularities of higher order appear at T = O. This behaviour is fully represented by equation (23). Equation (23) represents the real part of the Fourier transform of the one-sided PSD to' Oq.qs (eo); the corresponding imaginary part is its Hilbert transform according to equation (22). Once the cross-correlation function of the fictitious system subjected to a white noise input having a one-sided PSD has been constructed, then upon returning to the frequency domain according to the Wiener-Kinchine theorem, an equivalent expression of the one-sided WrOq,q,(w) which appears in equation (17) is obtained in the form 1 (j)roq.q,(w) = G',M'cW) =7T

foo R~~,( - T) e

iW T

-00

d r,

(25)

483

EVALUATION OF 8M

0'5

Of---\---f------+--/--H'--+-+-\-+~,.=-.~

-0,5

-10 -5 /. 2 r--------,----,------,.-----,- - - - - - , - - . ,

o

15 5 10 2r-----r---,---.,----..- -,.--,

(bl

-1·2

-2

5

6 (d)

(e)

if 2·5

4

0

0

-2'5

-4

-5

-15

-10

-5

a

5

10

f

W3 8 (T )

-6

-15

15

15

T

Figure 1. Cross-correlation function and its derivatives; Rarameters selected are (a) R~~~J,(T); (b) R~I~~,,(T); (c) R~I:~,(T); (d) R~~,,(T); (e) Rq~~,,(7").

Wk

== I, w, == 2, lk == ls == 0·1.

which can be rewritten in the explicit form

Upon recalling the main properties of the Fourier transform of the Hilbert transform of an arbitrary function b( T) [6,7], i.e., FT[HT{b(T)}]=-isgn (w)B(w),

(27)

where IT [ .] is the Fourier transform operator, sgn (w) is the signum function and B(w) is the Fourier transform of b( T), after some tedious algebra equation (26) can be written in an explicit form as

(28a)

484

M. DI PAOLA AND G. MUSCOLINO

1m [G2r,q, q, ( W )] = U(w)( -1)' [lHk(W )J2(E2r,skW + a2r,sk W 3) 4 -IHs(w)12(Ezr,ksw+a;r,ksW3)+ 2i

(_1)
(28b)

n=I,3,5

r>1

Re

[G2r+ 1,Qkq,( W ) ] = U(w)( _1)'+1 [IHk(W W(e2r+l,skW + arr+l,ks (

+ IHs(w )j2(E2r+l,ksW 2r-3

- L

(-1)

(n-l)/2

w2r -

n ,ks W

n]

a~r+l,sk(3)

3)

(28c)

,

n=I,3,5

r>1

1m [G- 2r+ 1,QkQ, ( W ) ] = U(w)( -1) r

[I

Hk(w)

1 2

('}'2r+1,sk W 2k

+ 8 2r+ 1,sk W 2 ) (28d)

The various parameters '}'j,b, '}'j,sk, Bj,ks> Bj,sk, Ej,l<.. and ej,sk (j = 2r, 2r+ 1) are given in the Appendix. Equations (28) show that the one-sided CPSD of the system subjected to a white noise input can be obtained as a linear combination of the squared modulus of each oscillator in decoupled form and in a recursive manner, whereas in equation (15) the mixed product of the transfer functions appears in a coupled form. It is worth noting that the real part of even CPSD and the imaginary part of odd CPSD have the same mathematical form and are the cosine Fourier transforms of the corresponding cross-correlation function. Moreover, a similar relation exists between the imaginary part of the even CPSD and the real part of the odd CPSD, which are the sine Fourier transforms of the corresponding cross-correlation function, Upon substituting equations (28) in equation (17) and integrating in the frequency range O-co, the SM in modal co-ordinates can easily be obtained in the form Re

_(_1)r v k v S[2 [A'2r,QkQ,] 2 AO,kk 'Y2r,sk W k

2r-4

+ A- 2,ss 8 2r,ks +):.

-

-

2

+ A2,kk 82r,sk + AO,5S'Y2r,ks w s

(-1)

nl2

] W2r-n-l,ksAn,ff,

(29a)

n-O,2,4 r>1

2r-4 + -A3,ssa2r,ks+ ):.

(-1)

(n-I)/2

]

W2r-n-l,ks An,ff '

(29b)

n-I,3,5

r>1

(29c)

485

EVALUATION OF SM

(29d) In these equations, An,if is the nth 8M of the input, given in equation (11), and Ar,J) (r = 0, 1, 2, 3; j = k, s) are the first four 8 M of the output evaluated by means of equation (12) and with the jth participation factor considered to be Vj = I: i.e.,

Ar.jj = f:
Ie

wrIH/wWOn(w)dw.

(30)

Once the 8M Ar,ii up to the third order have been evaluated, all other cross-SM of any order can be evaluated in a recursive manner by using equations (29). It follows that only 4m 8M instead of m(r+ 1)(m + 1)/2 8M must be evaluated by numerical integration, therefore drastically reducing the computational effort for the evaluation of all SM. Note that in equations (29) all the first terms are always convergent, while the last, involving the SM of the input An,if' can be divergent. As an example, if the input has a Tajimi-Kanai PSD (well known in stochastic seismic analysis) then the 8M of the input greater than AO,if are divergent. It follows that the 8M of the output greater than A4 ,Ql Q, are also divergent. In general, with exclusion of the case of white noise input, which has been studied in a previous work [4], if the nth 8M of the input is divergent the (n + 4)th 8M of the output is divergent too.

4. APPROXIMATE EVALUATION OF THE FIRST FOUR SM

As has been seen in the previous section, only the first four 8M having equal indices need be computed in order to evaluate all cross-8M of any order. In this section a numerical procedure is presented which, in a simple way, allows an approximate evaluation of the first few SM. To this purpose, in accordance with equations (11) and (30), the G.ff(w) function which appears in equation (11) and (30) can be written in the range O-w mx in the form N

0ff(w)= ~ Pk(W),

(31)

k~l

where N = w m x / Llw is the number of steps, W mx the appropriate cut-off frequency, Llw the frequency step selected such that the G.ff(w) is continuous in each step, and Pk(W) the interpolation polynomial defined in the form hew) =

(.~,-0 akiw') Wk(w),

(32)

where rk is the order of the polynomials selected (the maximum order should not exceed five [3]), ak are unknown coefficients, and Wdw) are window functions defined as

n ~ w:os; k- 1

elsewhere Here

nk = kLlw.

"l

(33)

Equation (32) can be rewritten as Pk(W) = nY(w )ak Wk(w),

(34)

486

M. OJ PAOLA AND G. MUSCOLINO

where n( w) and 8k are two vectors of order rk + 1. In the rest of this paper, the order of rk selected is three. In this case the vectors Ok and ak can be written as (35,36) The parameters can be obtained by imposing the continuity conditions on Gff(w) at the beginning and end of each discretized polynomial: i.e., (37) Here (38)

where the prime denotes the first derivative with respect to the radian frequency evaluated and ilk, the superscripts plus and minus can be avoided if the for the kth step in Gff(J2) is a differentiable one, and

at.,

(39)

By using equation (37), equation (34) can be written as h(w) = new)TNk1g k Wk(w),

(40)

It follows that the approximate PSD of the input can be obtained as Gff(w)=n(w)T[£ Nk'lgkWk(W)],

(41)

k=1

Once this has been evaluated, the SM of the jth single oscillator can be obtained by substituting equation (41) in equation (30): i.e., Ar,jj=

£ [fOO w !H(w )1n(w f W (w ) dW]Nk1gko 2

T

j

k=l

k

(42)

0

and upon using the definition of the window function Wk(w), equation (42) takes the form (43) where the vector k A ;'jj, which appears in equation (43), is given in the form (44)

where

s = r,. ", r+3.

(45)

Equation (43) shows that in order to evaluate the AT,jj, the 8M up to the (r+3)rd order must be computed, because the third order interpolation polynomial has been used. Moreover, as was shown in the previous section, the 8M up to the third order are sufficient

487

EVALUATION OF SM

to obtain any order SM of higher order. It follows that the various entries of the vector defined in equation (44) can be computed recursively as (see equations (29» A2

k -

..

r,JJ

=

(

-1

)

r [ k -

2

A O•1iY2r,jjWj

2r-4

1 + k A2,1i82r,jj+2 L (_ 1) 11/2 W2r-n-l •.ii k An

]

,

(46a)

n=0.2.4 r>l

(46b)

where k A"

n.

=

I

fl"+l_nn+l

wI! dw =

k

k-l.

(47)

n._, n+l The kAs'1i (s = 0,1,2,3) have already been evaluated [3] in a closed form solution, and for the sake of completeness are also reported here in the Appendix. Once the 8M up to the third order, defined in equation (42), have been evaluated, the modal SM of output can be evaluated by using equations (29), in which the SM of the input A",if are evaluated as follows: (48)

where (49)

Lastly, by using a piecewise polynomial interpolation of the PSD of the input, and using recursive formulas and the closed-form solution for the evaluation of the first few SM having equal indices subjected to the window function PSD, the approximate SM of the given system can easily be evaluated. It is worth noting that the only approximation made is in the representation of the PSD of the input, while the integrals present in equation (43) are given in a closed form solution. It follows that the sources of numerical errors are drastically reduced, and hence a very large step in the frequency domain can be used. 5. CONCLUSIONS

The spectral moments of MDOF linear systems, by using modal analysis, are expressed as a linear combination of the cross-spectral moments of the modal oscillators. The latter are the moments of the one-sided PSD of the input multiplied by the transfer functions of two modal oscillators. It has been shown that the cross-PSD of the various derivatives of the modal responses can be obtained in a decoupled form and in a recursive manner involving the PSD of each oscillator separately taken. In this way, the cross-spectral modal moments of any order can be computed from knowledge only of the spectral moments of each oscillator up to the third order. It follows that in order to describe the vibratory phenomenon, only 4m spectral moments (m being the number of the first few modes selected for the analysis) need to be computed by numerical integration. The latter quantities have been computed by using the method suggested by Pfaffinger [3], improved by the use of the recursive relationships presented here. REFERENCES 1. E. H. VANMARCKE 1972 Journal of the Engineering Mechanics Division, ASCE 98, 425-446.

Properties of spectral moments with applications to random vibration.

488

M. DJ PAOLA AND G. MUSCOLlNO

2. A. DER KIUREGHJAN 1980 Journal of Engineering Mechanics Division, ASCE 106, 1195-1213. Structural response to stationary excitation. 3. .0. O. PFAFFJNGER 1980 Computer Methods in Applied Mechanics and Engineering 24, 269-286. Analytic evaluation of modal covariance matrices. 4. M. OI PAOLA and G. MUSCOLINO 1986 Journal of Sound and Vibration IUl, 233-245. On the convergent parts of high spectral moments of stationary structural response. 5. D. O. PFAFFINGER 1983 Journal of the Engineering Mechanics Division, ASCE 109, 357-372. Calculation of power spectra from response spectra. 6. A. PAPOU LIS 1965 Probability, Random Variables and Stochastic Processes. Tokyo: McGraw-Hili Kogakusha. 7. M. OI PAOLA 1985 Solid Mechanics Archives 10,225-243. Transient spectral moments of linear systems under random agencies. APPENDIX

Coefficients in equations (28) and (29). The various parameters which appear in equations (28) and (29) are given as follows: 'Yj.•k = (-OJ[aj,'skCkwk+.a ),skWk],

'Yj,k. = a tu», + .a!.ksWso

(Ala, b)

8j ,sk = (-1)l[ aj,'.dkwk -.aJ.skWk],

OJ,b =

a ;k.C.W s - .a!.ksW.,

(A2a, b)

Ej,sk = (-l)j[2CkWkWk.aJ:sk+ aJ~sk(w~ -2w~)],

(A3a)

Ej.ks = 2.alks!:,w,w s + a;k.,(W~ - 2w;).

(A3b) These parameters are evaluated by computing the coefficients of the Fourier transform of the correlation function defined in equation (26). First few SM having equal indices. Here the first few SM having equal indices subjected to real random processes having a PSD window function of the type defined in equation (33) are given in a closed form solution. These quantities have been obtained analytically by means of partial fractions [3,5]. As has been shown in the previous section, upon knowing kA r .1!- defined in equation (44)-and evaluating it for r = 0, all other kAr,jj vectors can easily be obtained in a recursive manner, To evaluate kAO,jj, the following relationship can be used:

.» = X/d j •

(A4)

k Ao

k

dj is a vector whose entries are given as kdlj =! In (kb 1) + (Wj!CJwj){tan- 1 (kb2j ) + (7T /2)[1- sgn (kb2j )]),

kd2j = (li!:jCJ)j) {tan -1 (kb2j) + (7T/2)[1 - sgn (kb 2j)]), kd4j = 0/ !:jWj) tan- 1 (kb4j), kd3j =1 In (kb3j) - (Wj/ !:jWj) tan- 1 (kb4j),

(A5)

where kb _ ({}k- Wj)2+CJWJ lj-(n -)2 +~jWj y2 2' Hk-l-Wj kb _ ({}k+ Wj)2+C]W] 3 r (fl k- 1 + Wj)2 + !:]CJ)J'

(A6)

and the matrix Xj is given as X.= 1

- 1/ 4WJWj 0 1/4Wj [

!

1/2w] 1/4wj 0

1/4wfWj 0 -1/4wj

-w]/4wj

1

1/2WJ] -1/4Wj 0

.

(A7)

w]/4wj

Equations (AS) are quite different from those reported in references [3,5] because the presence of the term (7T /2)[1- sgn (kb2j)], which is essential for completeness.