On the response of single and multidegree of freedom systems to non-stationary random excitations

On the response of single and multidegree of freedom systems to non-stationary random excitations


2MB Sizes 0 Downloads 9 Views

J. Sound Vib. (1968) 7 (3), 393-416



J. K. HAMMOND Institute of Sound and Vibration Research, the University of Southampton, Highfield, Southampton, SO9 5NH, England (Received 26 September 1967)

The early part of the paper indicates two current theoretical spectral methods of dealing with non-stationary processes. Then, by using Priestley’s evolutionary spectral density, the non-stationary response of single and multimode systems to stationary excitations is derived. Finally, by constructing a non-stationary evolutionary spectral function, the corresponding non-stationary response of linear systems is investigated.

1. INTRODUCTION The study of stationary random processes within an engineering framework has been covered exhaustively and the formal methodology for dealing with such processes is now well known. The mathematical descriptions of the physical phenomena involved can be studied either in the time or frequency domains. These are equivalent and are related by the Fourier transformation. The choice of domain is merely one of mathematical convenience. The temporal description of the process can be regarded as the “natural” description whilst the frequency description is a measure of the process power distribution over frequency, i.e. the process is assumed to be composed of a sum of a continuum of sine and cosine waves of different frequencies having a random amplitude distribution. For a process x(t), this is expressed as x(t) = ?’ efwf dX(w) -co


where X(o) is a distribution function and is orthogonal, i.e. random variables dX(o,) and dX(wJ are uncorrelated when w1 # w2. This leads to the power spectrum defined as an ensemble average of the power at each frequency: dF(o) = E{jdX(o)j’j, where E denotes the mathematical expectation. It is then useful to conceive of a power spectral density functionf(w)

(2) where

This spectral description is very useful in engineering since it has the direct physical interpretation of being the process energy density distribution over frequency. Also, when ergodicity may be assumed, spectral density functions can be obtained from measurable 393



quantities. The direct spectral functions can be estimated either directly by using narrow band filters, or by Fourier transforming the autocorrelation function. Cross-spectra are obtained by taking the weighted Fourier transform of the appropriate cross-correlation functions. Further, this description fits neatly into place in the analysis of the response of linear time-invariant systems to stationary random excitations. For example, the mean square response of a linear system to a stationary excitation is an infinite integral over frequency of the output spectral density, which for a single-degree-offreedom system is given by f,(m) = Ia(w)12fx(w),


where &(w) is the input power spectral density, f,(w) is the output power spectral density, and a(o) is the system receptance function. Multimode systems have more elaborate but fundamentally the same form of description. However, it is an unfortunate fact that most random processes in nature are non-stationary. An immediate consequence of this is that ergodicity must be rejected. This means that the statistical descriptions normally employed are time dependent, and at once the computation of averages of any sort is a difficult task involving a large number of records. However, by making some assumptions about the process, meaningful analyses can be conducted which do not involve prohibitively large numbers of calculations. The functions used to describe non-stationary processes are generally extensions of familiar stationary functions but it is their interpretation that is important. In many cases the number of records available from a process may not be large and in fact there may only be a single record (e.g. missile launch, where cost per launch is high). Approaches have been devised to try to deal with these cases as well. This implies that in all subsequent records the non-stationarity will have the same form relative to some starting time. This is an assumption of an ergodic nature, i.e. that an average or averages over one record are representative of averages taken over hosts of records at a fixed time instant. This is a reasonable assumption in a large number of physical situations. As noted earlier, dynamic system response is the main application of random data studies and so, after a brief description of a few approaches to the problem the response of linear time-invariant systems is considered within the framework of one of the more recently proposed techniques. 2. GENERALIZED SPECTRA The generalized or double-frequency spectral density function is a direct extension of the function used in stationary process analysis. It is a mathematical procedure and has little physical significance. Its proper description is that it is a double Fourier transform of the correlation function. The non-stationary cross-correlation function is defined as

and the corresponding generalized spectral density is given by m


w2) = ss


R,,(t,, t2)ei(“l’l-“*t2)dfl dt,.


The inverse relationship is &,,(t,, t2) =

f&,(co,, w2) e-i(W~‘~-WZfZ’dw, dW2.






This spectrum can be discussed a little further by considering spectral representations the form x(t) = f eiWtdX(w). -CO


This type of process is termed “harmonizable”. Here X(W) is no longer an orthogonal process since x(t) is non-stationary For a(t) real we can rewrite (8) as

[l, 31.

x(t) = f eeiUf dX*(w). -_g) Here * denotes the complex conjugate. Combining equations (8) and (9) to form &(t,, R&t,, tz) = E[x(t,)x(t,)]



t2) we have

= f JPE[dX(w,)dX*(w2)]e-i(Y1’1-~zr2). -cc --Lo


Comparing equations (10) with (7) we see that fxx(~,,~~)du~

do, =



This is a correlation in the frequency domain of the amplitude functions, which in the stationary case are orthogonal processes. 2.1. INPUT-OUTPUT





Although a direct physical interpretation of the double-frequency spectral function does not exist it can be quite useful analytically in obtaining the response of linear, time-invariant systems. The non-stationary input-output relationship, developed similarly to the stationary case, becomes, by virtue of the “non-orthogonality” property [l, 21, .&&JI, ~2) = a*(@,) a(wz).L(w~, ~2)


for a single-degree-of-freedom system. Here fY,,(w,,wZ) is the output generalized spectral density function, f,,( w1,w2) is the input generalized spectral density function, and a(u) is the system receptance. J. B. Roberts [2] has used the generalized spectral density to deal with the class of processes termed evolutionary white noise : those with correlation functions of the form K&l> t2) = WI) WI - f2) where I(tI) contains the time dependence. Roberts has shown that the response of any system to a periodic I(t) may be determined graphically from a vector plot of the system’s receptance function, which could be obtained experimentally. Further interesting results involving generalized cross-spectral density functions are obtained in reference 1. Their direct usefulness is again subject to the ease of calculation of the relevant spectra. 3. EVOLUTIONARY


M. B. Priestley [3,4] has defined a non-stationary spectral function which is closely related to the conventional spectral density, so as to preserve a physically meaningful interpretation.



He called it the evolutionary spectral density function, and showed that, under certain conditions, it can be estimated from a single trace. Priestley explains both in rigorous mathematical terms [4], and then more intuitively [3], the underlying principles justifying, and also serving to physically interpret, his choice of spectral representation. What follows in this section is a brief, non-mathematical summary of what is elegantly put forward in references 3 and 4. Both stationary and non-stationary processes are representable by the integral equation

~0) = / 40, w> dZ(w)



where Z(w) is an orthogonal process and $(t,w) some suitable function. For a stationary x(t) one possible form of 4(t, CO)is eiw’. When x(t) is real, the conjugate of (13) is x(t) = r $*(t, w) dZ*(w). --m


Multiplying (13) and (14) we get


r ~C(t,~~)C*(s,W2)dZ(W,)dZ*(W (1%



Taking expected values and putting $(t, CO)= eiw’for a stationary process, then

-w~)x(~)l= Jp 1 eiwl’e-‘“2sE[dZ(o,)dZ*(o,)]. --m--m

Since E[dZ(w,)dZ*(w,)]

= S(w, - w,)EldZ(w,)l*

for t = s this becomes Varx(t) = E[x*(t)] = f E]dZ(w)]*. --m


Varx(t) is a measure of the total power of the process at any instant and can be determined in terms of a power spectrum as Var x(t) = f d&‘(w), --m


dF(w) = E]dZ(w)l*.



Now, assumingf(w)

exists, we write: dP(w) =f(w) dw,


where f(w) is the power spectral density function. So the choice +(t,w) = eiWt yields the decomposition of the process into sine and cosine waves of different frequencies, and the spectral function so obtained can be physically interpreted as an energy distribution over frequency. Priestley notes that it has been shown that the existence of a particular #t,w) means that there exists an infinite number of functions that allow representation (13). Therefore some form other than &t,w) = eiwt could have been chosen for the stationary




case, with its corresponding orthogonal Z(o). However, the new spectral density could not be interpreted as an energy distribution over frequency (in this case w being merely any variable). When dealing with non-stationary processes the choice (b(t,~) = eiWtis inadmissible for a corresponding orthogonal process Z(o). Relationships can be developed using I$(&w) = eiwf but due to the resulting non-orthogonal character ofZ(w), two-dimensional spectral functions must be considered and little physical significance can be attached to these (see section 2, above). It is clear therefore that some other form of +(t,w) is desired, and to preserve the notion of frequency in the resulting spectral density Priestley introduced his so-called “oscillatory” function +(t, ~0)= A(t,W)e*B(w)t which can be looked on as amplitudemodulated sine and cosine waves of frequency O(U). This function preserves the idea of frequency as long as the modulating function varies slowly with time, i.e. the modulus of the Fourier transform of A(t,w) with respect to time has an absolute maximum at zero frequency for any w, i.e. A(t,o)

= j eitu dH(u, o), -al

IdH(u,w)] having an absolute maximum at u = 0. Assuming o to be single valued so that the variable w in (13) can be transformed to 6(w), then by redefinition of Z(w) and A(t,w) we can write x(t) = f A(t, w) [email protected]”dZ(w). --m


The case A(t, UJ)= 1 and 0(w) = w (stationary process) is included in the class of oscillatory functions. As before, using the orthogonality property of Z(w), we have E[x2(t)] = Varx(t) = f ]A(t,w)j*dF(w), --m


where dF(w) = EjdZ(w)j*. Varx(t) is a measure of the total energy of the process at time t. We can write (21) as Var x(t) = 1 dF(t, w), --m


where dF(t, w) = jA(t, w)]*dF(o),


and is seen from (23) to describe the process energy distribution over frequency. Now defining density functions dF(t, w) =f(t,

w) dw


dF(o) =f(w) dw,


we see that f(t, w) = IA& @)I‘f(w).


This is the evolutionary spectral density at time t and represents the energy spectral content of the process near t. Priestley [4] discussed at some length the characteristics of the modulating function



A(t,o). As there is an infinite number of possible choices of c$(~,w) [and therefore A(t,w)] x(t) is representable in terms of functions which have different rates of change. Since, as is shown in reference 4, the estimation of the spectral density function f(t,~) involves an average over time it would be best to use the particular A(t, w) whose temporal fluctuations are slowest in order to get the best possible estimate of f(t,~). Priestley shows in his estimating procedure that it is not necessary to specify any 4(&w), but that the estimation procedure involved automatically selects the optimum A(t, UJ).Therefore explicit knowledge of A(t, w) is unnecessary in the estimation of an evolutionary spectral density.

4. AN OSCILLATORY PROCESS In practice there are many situations where non-stationarity is significant and the question arises as to which processes are amenable to Priestley’s analysis. As one example of a non-stationary process we shall examine the case of a pressure field close to a jet engine. We shall consider the changing spectral patterns obtained near an engine as the jet velocity varies in the vicinity of some arbitrarily chosen figure of V,, ft/sec and attempt to interpret

Figure 1. Measurement position P relative to jet. the resulting varying spectral distribution as an evolutionary spectral density function. We shall not be considering the detailed mechanism of sound production here; indeed the discussion that follows is intuitive and empirical and is based on experimental results. The experimental results used are soon to appear as Royal Aeronautical Society Data sheets under the title “The estimation of near field sound pressure levels due to jet noise.” We shall use the data sheets to calculate the acoustic spectral levels due to a typical jet engine running with mean fully expanded jet velocity Vj ft/sec, at a point P, fixed relative to the exhaust (Figure 1). V, will take a range of values in the vicinity of V,. The resulting spectra will, of course, be the stationary spectra for the particular conditions. We shall choose typical numerical values, e.g. V,, = 2000 ft/sec, a jet nozzle diameter D = 2 ft, and a static jet pipe temperature T,, = 900°K. Vj is chosen to vary from 1700 ft/sec to 2000 ft/sec (=V,) in steps of 50 ft/sec. The seven spectra obtained at position P (Figure 1) are calculated for a frequency range of 100 Hz to 1000 Hz. This encompasses typical aircraft panel frequencies. Now consider the practical case of the jet running in the vicinity of 2000 ft/sec. For example, the throttle is closed from 2000 ft/sec to 1700 ft/sec and then opened again within a certain time duration. A single pressure transducer measurement is made at P during this procedure giving a single pressure-time trace. Applying Priestley’s method of analysis to the record we would hope to attain, if the process were oscillatory, approximations to the stationary spectra already obtained above, at the appropriate times (velocities). The





question is whether the process is oscillatory. We are assuming here that the time varying conditions in no way alter the fundamental physics of the situation. Priestley points out that an oscillatory process is representable by x(t) = i A(t, w) eiwf dZ(w), --m with Z(w) an orthogonal process and A(t, w) having a Fourier transform that has an absolute maximum at zero frequency. This leads to the evolutionary spectral density f(C w) = [email protected], w) I*fb), where f(w) = dF(w)/do. He also states that there are a large number of different functions A(t, w) with corresponding differentf(w). Now let us choosef(w) here as the stationary spectral density obtained from the jet running at 2000 ft/sec [i.e., f(w) =_&,(w)]. Being stationary this is composed of an orthogonal Z(w). Now if we can find an IA(t,u)12 such that IA(t,w)12*fV,(w) gives good agreement with the stationary spectra already obtained at the relevant times (velocities), and also if A(t, w) satisfies the previously stated conditions, then the process is oscillatory. The determination of ]A(t,~)l~ involved a curve fitting procedure that was aided considerably by the form in which correction terms are included in the data sheets. This yielded IA(t,w)l’=(l




u(t))2 22222.2

where V0 is the “datum” jet velocity (2000 ft/sec), u(t) contains the time dependence and is the decrease in jet velocity giving the true jet velocity as (I’, - v(t)), f is the frequency (=w/27~), andf, is a constant frequency = 365 Hz and is apparently a function of the position P relative to the jet. I












11~51 0









200 300 400 500 600 700 Frequency




800 900 1000



Figure 2. Jet noise spectra for u(t) = 0 and v(t) = 50 ft/sec. The continuous curves are from the data sheets. The triangles show the values obtained by multiplying the curve for u(t) = 0 by the factor ][email protected],w)] 2 with v(t) = 50 ft/sec. Sample results obtained using this expression are shown in Figures 2 and 3 together with the curves obtained treating the process as stationary. The agreement is seen to be good. It is also noted that the process is not uniformly modulated, i.e. lower frequency spectral levels are augmented, whilst higher frequency levels are scaled down for lower jet velocities, The “pivot” frequency is 365 Hz. This effect is small, however, in comparison with the overall lowering of spectral level with decreasing velocity. We choose A(t,w) to be the real 21



positive square root of equation (26) and if this A(t, w) satisfies the required conditions then the process is oscillatory. The velocity variation with time will be chosen as a linear variation for simplicity, symmetric about t = 0, i.e. v(t) = kt t>o, = -kt I




t (0. I






12,0I I.9 .!


E p







I 100

I 200

I 300

I 400

i 500













Figure 3. Jet noise spectra for u(t) = 0 and u(t) = 250 ft/sec. The continuous curves are from the data sheets. The triangles show the values obtained by multiplying the curve for v(t) = 0 by the factor [A(?,w)12with u(t) = 250 ft/sec. This means that the factor (1 - v(t)/V,)* will = 0 for f > V,/k and for i G -V,/k, that the process only exists for -VJk G t G V,/k. Therefore


which is real valued, non-negative and being zero for 1t 1 > V,/k has a Fourier transform. Thus the process is oscillatory. By specifying f(w) and hence the particular orthogonal Z(W) we are probably overrestricting A(t,w) and there might be a large number of functions with wider “characteristic width” [3,4] than the one we have obtained. The main physical feature of the jet noise problem is the velocity of the jet relative to the surrounding air, i.e. (V, - v(t)). The problem could therefore be looked on as being representative of the case of an aircraft rolling along the ground with constant exhaust velocity of V, ft/sec and varying ground speed. The temperature correction term is now unnecessary and the expression for A(t, co) would be something like 4&%~()

-(l-f&J(C) 5300x103



It may not be quite the same as (28) owing to the change in the physical situation but is certainly likely to be very similar. Therefore, included in the processes for which time-varying spectra can be computed by Priestley’s method are likely to be noise analysis for most aircraft and missile test and flying conditions. The particular A(t,w) obtained here is used in section 7 below in response calculations.







An engineer’s interest in statistical descriptions of random data is centred on their application to commonly occurring problems in order to get better estimates of the effect of random processes. We will now consider the response of linear, time-invariant vibratory systems to random excitations. Having a formal approach to the description of a class of non-stationary processes due to Priestley, we wish to investigate how it fits into the framework of conventional frequency domain response work. The evolutionary spectral density has a close relationship to the familiar stationary spectral function and so it is likely that the response equations will be similar. They will, of course, be functions of both frequency and time. Before considering the problem of the response of a linear system to a non-stationary excitation we will first consider a far simpler problem. This is the problem of non-stationary response of a linear system to a stationary excitation in which the non-stationarity in the output is due to the build up of the response from rest. This problem has been analysed already using a different technique by T. K. Caughey and H. J. Stumpf [5,6] for a singledegree-of-freedom system and later by Y. K. Lin [7] for multimode systems, so solutions are available for comparison. We will first consider the response of a single-degree-of-freedom system. 5.2.


Let us examine the response of a second order, linear, time-invariant system to a stationary excitation commencing at t = 0, i.e. the excitation for t za0 is a portion of a stationary process. For a time the response will be non-stationary as the response builds up from rest, finally settling down to be stationary. This is on the assumption that the system is asymptotically stable. Let 5 denote the damping ratio and o. the undamped natural frequency. Then the equation of motion is ji + 25w,j + w$v = x(t), (29) where x(t) = 0 for t -c 0 and is a stationary random variable for t > 0. Since x(t) is stationary it can be written x(t) = jp eiwrdXO(o).


--co The superscript “0” denotes a stationary process. X0(w) is orthogonal. The response y(t) is non-stationary and let us assume that it is representable by y(t) = 7 G(t, co) eiWtd Y(w),

t > 0.


--m Y(w) is an orthogonal process and G(t,w) is the modulating function making y(t) an oscillatory process. Substituting (31) and (30) into (29) and taking the time derivative under the integral sign we finally obtain [writing G for G(t, w)], m {G + 2G(iw + two) + G(wi - w* + 2i&owo)} eiWtd Y(o) = f eiW’ dXO(w).

s --m





Now consider the equation (29) for large values of t. For 5 > 0, y(t) will ultimately be a stationary random variable effectively representable by m

y(t) =

1.eiwt d YO(,),



Substituting (33) and (30) into (29) we get x(t) =


= f (CIJ~ - w2 + 2i{woo) eiWtd Y’(w).

We shall write (wi - w2 + 2icwoo) as l/cc(w). Now forming the expression for x(t,) as


x(t,) = f e-‘“‘ldXO(w)* = 1 m


e-iW*ld YO(,)*,



and multiplying (34) and (35) and taking expected values, we obtain the autocorrelation function of x(t) as &(G tr) = jp elCwtPwltl)E[dXO(w) dXO(w,)*] --m m



c&J) c+,)*

ei(wt-wlrl)E[d YO(o) d Y”(wl)*].



Noting that function X0(w) is orthogonal we write dF,O(w)= E (dXO(w) 12. Then assuming the existence of a spectral density function we write dF,O(w)=fz(W)dw. Doing the same for Ye(w) permits equation (36) to be written as m








and the two descriptions of the autocorrelation (t - t,) only. By the uniqueness of representation

function of x are seen to be a function of of a function as a Fourier integral we can

equate integrands of (37) to get Ia(~)l’fX~)



This is the well known result for a single-degree-of-freedom system. We could have obtained (38) alternatively by equating squared integrands in (34), i.e. we equate the power at each frequency. This gives ]dX”(w)j2 = ‘1 ‘(“(ii”. aw

On taking expected values, since we must average the power distribution over all records to obtain a meaningful description of the process as a whole, we obtain (38) once again. We shall attempt to deal with equation (32) similarly, remembering that the evolutionary spectral density for a process given by Z(t) = f A(t, w) eiWtdZ(w)






isf=(t, w), where fi(&w) = IA(0J)12fi(W).


We shall rewrite (32) for brevity as eiwt dX”(w) = x(t)


--m where TV= G + 2G(iw +
2ldY(w)12 = IdX”(o)12. Taking expected values and assuming the existence of spectral density functions we obtain (43) But by (40) the evolutionary spectral density of v(t) multiplying both sides of (43) by I G(t, W)I2 we get


f,(t’w)= Ip+


is fy(f,w) = IG(t,w)12&,(~),



This can be rewritten as (45) This is an expression for the evolutionary the system receptance and a time-varying contains derivatives of G and as t --f DJthe these derivatives -+ 0, so that for large t we

output spectrum in terms of the input spectrum, portion due to system transients. The term p system tends to a “steady-state” condition and can write (45) as

which is as we would expect. However, we are interested in what happens for t near zero and so would like to be able to evaluate the factor jpt~/G + 1Im2exactly at any time t. This will be attempted by using the argument discussed below. We know, by equation (40), that .A(&w) = IW, w>I 2fv(w>, where f,(w) is obtained from an orthogonal


Y(w). Ultimately, however, (38) will hold since

y(t) will be stationary. Let us choose the f,(w) of (46) to be f ,"(u). Now we must find the corresponding G(t, co). Since we have so restricted f,(w) and also since we know that

AJ~Q9:,“) -+fvohJ) then from (46) we see that G(t,w) must obey IG(t,w)l’+





Also the equation (31) is restricted to using d Y(w) = d Y”(w), and so from this equation we see that G(t,o) -+ 1. f--Ku

Thus we know the asymptotic behavior of G. Since we know the ultimate (steady-state) value of G then the coefficient of eiWtdY(w) in (32) [where the d Y(w) is now d Y”(w)] can be looked on as the left-hand side of a differential equation in G, the right-hand side of which must be arranged to give G = 1 for large values oft. This differential equation is G + 2G(iw + {wO)+ G(w; - w2 + 2i5ww,,) = (w; - w2 + ~z&Jw,).


We shall use (47) to rewrite equation (45). Equation (47), in the notation previously used, is

and therefore I-G



&JJ) Substituting (48) into (45) we get f,(t,o)


= I~(w)12L%) IG12.

So in order to obtainf,(t,o) we must solve (47). Since the system is quiescent for t < 0 the initial conditions we use are G(0, w) = 0, G(0,0) = 0. Using solution G(t, w) = A eht and the conditions above we obtain G(t,w)=


co0 + i(w + &) e-l;wote-i(w-8)

t +


t e-i(w+S)




where w = wod(l - 52). We now want IG12= G. G*. This comes out to be po;+ClJ2-cijz 1 + e-2~war[1 + 2


- 2e-~wot<~cosWtsin&t w

- 2e-+“‘coswtcos8t


1 -


- 2 Ff e-cwot sin wt sin 63 . w 1


The mean square of y(t) is therefore

Eb2(t)l = ~_Uf,w)dw = ~f,o(w)la(w)12lG12dw




where IG12is given by equation (51). This result [equation (52)] agrees identically with the result obtained in (5) and (6) where the mean square response is obtained from the expression [using the notation of (6)] f--f0 t--fly

HY~(~)I= a:(t) = K& t>= / dt, s d52WA%) k,(t - 51,t 0







k,,,,(f,, t2) and k,,(t - [i, t - f2) are covariance functions and h the impulse response function.

When x(t) is stationary

the stationary auto-correlation function. Rx, is replaced as the Fourier transform of the stationary spectral density function and then the two temporal integrals are evaluated for a finite operating time for a second-order system to obtain the final result for a:(t) in terms of an integral over frequency. We have achieved the same result preserving the notion of the evolutionary spectral density throughout. 5.3.







We shall now extend the analysis for a single-degree-of-freedom system to include non-stationary multimode system response to a stationary excitation. The normal mode approach is used as proposed for example by Powell [8] for stationary excitations and responses. Let the equation of motion for a distributed parameter system be mti +

C+ +

Q(w)= p(r, t)


where Q is a linear differential operator in space and r denotes the spatial position at which the pressure acts. If a solution is assumed of the form

where &(r) is the mode shape in the nth mode and qn(t) the generalized co-ordinate corresponding to the nth mode, then it can be shown that the equation of motion can be written, for motion in the nth mode, as (54) where M,, C, and K,, are the generalized mass, damping and stiffness, respectively, and are defined by M. = s &Jr)



C, = s cf’,(r) dr A

and K, = CIJ~ M,, where w, is the natural frequency of the nth mode. L,(t) is the generalized force in mode n and is defined as Ut)

= 1 ~(r, t)f.(r) dr.



We can rewrite (54) in the form MA& + 25, w, 4” + w; q,,) = &,(t) where 5, is the damping ratio of mode n.




Once again, the system is assumed quiescent till t = 0 when the motion is allowed to build up in response to a stationary excitation p(r, t), For an asymptotically stable system and for large values of t w(r, t) [and so q,,(t)] will be stationary, so we can write

J*[email protected])PdQOb) --m

w(r, t) = i eiw’dWo(r,w)=ng



where the superscript “0” denotes stationary processes. Substitution of eiwt dQ,O(w)

--m into (56) and use of (55) gives m M,{-w2 + 2i(&ww, + wt } e '"'dQ~(~)=lp(r,,t)f.(r,)dr,. (58) s --m A The conjugate equation of equation (58) in mode s, writing Z,(w) for MS{-w* + 2i{,




is m J Z,(w) * e+“” de:(u)* --oo

= /P(r,,

t,)fs(rJ dr,.



Multiplying (58) and (59) and taking expected values we obtain


J pZ,,(~)Zh)


-a, --m


s.iptr,,t)P(r2,~l)_Ltrl)fstr2)dr, dr2.

If we write eiwt dP”(r, w), --m

the above equation becomes meiwlt-tl)Z”~~)Z~(,)* E[ dQi(w)dQt(o)*] I -03 m

{E [dPO(r,, co) dPO(r,, co)*] eiw(‘-“)fn(r,)fs(r,)> =SSI -co

dr, dr,.

A, A2

Let and let

exist. This is interpreted as a modal cross-spectral density. Similarly, E[dP”(r,,w)dPo(r2,W)*]





dF”(rl, r2, w) dw ’





is the spatial cross-spectral density of the excitation (which for a homogeneous field is a function of jr, - rzj only). Therefore equation (60) becomes .? J -m


= r 1 1 e’w”-tl’f,(r,)fs(r,)~-~(r,,r,, -co A, AZ

co)dr, dr, dw,

from which, by the argument previously put forward, we can write


or w(r, t) = J” eiWtd WO(r,w) = ngl f,(r) 7 eiwt dQt(w). --m --co


Equating squared integrands here we get IdWO(r,w)l’=

5 n=l

5 Sn(r)fs(r){dQnO(W)dQsO*(W)}. s=l

On taking expected values and assuming the existence of density functions,

Substituting for f:(w),, from (61) into (63) we obtain

fn(rlMr2)[email protected]~r2, 4 dr, dr2


AI AZ wheref$(r,w) is the displacement spectral density function at r. This is a simple form of the result obtained by Powell for a stationary excitation. Now we shall attempt to derive the non-stationary response equations due to transients such that the asymptotic behaviour of the general time dependent equation is identical to equation (64). In this case since the response is non-stationary let w(r,t)=

r B(t,w)e’“fdW(r,w)



q&) = 1 GbW”‘dQk4


so w(r, t) = J B(t, w) e*Wtd Wr, w> = ng, f,(r) J’ G(t, w>eiwrdQ,.


Equating the power at each frequency






Taking expected values and assuming the existence of density functions we get f&r,

a) = g


5 L(r)_&(r) G,(t, w) G%, w)_L&J)~~.



Now we shall restrict ourselves to using random processes whose ensemble average gives i.e. that spectral function obtained in the steady state. This means we restrict fQoW"S, G,(t, w) to become unity for large values of t. We already have the expression for fz(o)nsin equation (61). So substituting (61) into (66) we obtain the following result for the evolutionary spectral density of the displacement at point r :

In order to solve (67) it remains to evaluate G explicitly. We recall that q,,(t) = i G,(t, o)efWt dQjl(w). -cc Substitution of this q,,(t) into (56) gives mM.{Gn+ 2G&w + ~,,uJ,,)+ G,(wi- t.2 + 2i5,wwn))eiwfdQH(~>

f --m


= m




As before, we impose the condition that the asymptotic value of G, is unity so the quantity within parentheses in (68) can be looked on as the left-hand side of a differential equation in G,,, the right-hand side of which must be arranged to give G, = 1. t+m

This equation is G” + 26&w + L,w,J + G,(w: - CJ + 2i1&4J

= (of - w2 + 2i[,ow3.


Using the initial conditions G,(O,o) = 0, G;,(O,w) = 0 since the system is quiescent till t = 0, we solve the equation and substitution of the solution into (67) gives the final result : co

f,(t, r, ~1=



cc “=,



_fXrl_fXr) =y

fn(rllfs(rdf%-,, r2, 4&dr2


)I (cos8,t+(i”W1+iW)sin8,f , 0s 11 wn



1 - e-cswst eiwt coscii, t +



(5, w, - iw> sin&t _

Gil,= w,(l - z;i)““. The mean square of w(r, t) is Hw2(r, t>l where f,(t, r, CLJ) is given by (70).

= 1 _A&r, w>do --m






This result agrees identically with the result obtained by Lin [7]. It is in fact a particular case of Lin’s result since he calculated the non-stationary cross-correlation function of the response E{w(r,, t,)w(r2,tZ)} whereas we have considered the mean square response at a point, i.e. r, = rz = r and t, = t, = t. As with Caughey and Stumpf’s treatment of the single-degree-of-freedom system Lin obtained his result via a finite time integral of the Green’s function for the structure whilst in the analysis the idea of a time-dependent spectral density has been preserved throughout. 6. NON-STATIONARY RESPONSE OF A LINEAR, TIME-INVARIANT SINGLEDEGREE-OF-FREEDOM SYSTEM TO A NON-STATIONARY EXCITATION We shall now consider the non-stationary response of a linear, time-invariant singledegree-of-freedom system to a non-stationary excitation. The non-stationary regime of the excitation will be assumed to be embedded in an otherwise stationary process, i.e. the record of the excitation may be considered representative of a stationary process for all times except for t lying in some specified interval where the process is known to be nonstationary. The process will be assumed oscillatory. The response equations will be developed using arguments similar to those put forward in section 5 above. The preceding rederivation of results already available in references 5 and 7 is a source of some confidence in applying the same method to the problem involving a non-stationary excitation. The equation of motion we shall consider is

j; + 25wi3j + w;y = x(t)


where x(t) is now an oscillatory process so x(t) = 1 A(t, w) eiwt dX(w), --m and as before y(t) =

1G(i, co)erwr d Y(w). --m

Substituting these equations into (71) gives x(t) = J (G + 2G(’zw + &.o,,)+ G(wi - w* + 2i5wwo)) eiWrd Y(o) --m -


A(t, w)efWr dX(w).


s -03

This gives two representations for x(t) as the limiting force of the sum of sine and cosine waves with different frequencies and time-varying random amplitudes. For convenience we Put D(t, co) = [G + 2G(iw +
1A(t, w) ei”‘dX(o) --co

We shall consider two representations equating squared amplitudes,

= J’ D(t, w) efw’d Y(W). -CO whose power at each frequency is the same, i.e.

I&t, w)l* ldX(o)l* = 1% o)l* IdX(w)l*.



If we take expected values and assume the existence of density functions, remembering that D s p + G/a(o), then (73) Now the evolutionary by [Gj2 we see that

output spectral density is fy(t,~) = IG]‘f,(w), so multiplying (73)

or fy(C w) =X&9 ~)l&J)l”.

jpa(U);G + li2 9


wheref,(t,w) is the evolutionary excitation spectral density. The factor I~LQ(o)/G+ 1le2 contains the contribution of the system’s time-dependent response. This must be evaluated to findfv(t, CO)and we shall now attempt to do this. Priestley [3] makes the point that for a particular process the oscillatory function in the integrand can take different forms. However for convenience here we choose our function D(t,m) to be a function that varies in exactly the same way with time as does A(t,w), i.e. we say that g;

= h(w)

where h(o) is some function of frequency. Now in order to determine h(w) we restrict our non-stationary effectively stationary for large times, i.e.


excitation to one which is

A(t,w)+l. *‘CC

This means E 1dX(w)] 2 =fz(o) d w where f:(w) is the stationary spectral density of the process at large times. Since the system is positively damped, G finally becomes a constant which we again choose as unity so that E Id Y(co)/~=f,O(w)dw wherefz(w) is the stationary output spectral density. Since A(t,w) and G(t,w) are effectively unity for large values of t then D(t, w) -+ -!t-%0 da) so h(w) = l/a(w).

Thus we have the relationship D(t,o) -=-* [email protected],0)

1 +.J)

Writing this out fully we have G + 2G(io + SW,,)+ G(w; - w2 + 2i5wwJ = A(t, co) (co; - w2 + 2i5wwJ.


So by assigning a particular form for D(t,co) and defining the asymptotic behaviour of conveniently we have obtained a differential equation for the modulating function G for process y(t). This equation reduces to (47), for the case when A(t,w) = 1 for all t (stationary excitation). A(t,w)





Substitution of (76) into (74) gives fy(GW) =r,(r,w,l+)12~~2.



So to obtainf,(t,w) we must solve (76). We shall now develop the analogous equations for a multimode system. 6.1.








We will examine the response of a multimode system to a non-stationary, spatially distributed excitation. The non-stationary time trend is assumed to be the same at all spatial positions. The arguments used are similar to the previous section. Now

Let u’(r, t) = J’ B(t, w) eiwt d W(r, w), --m

q&l =


1GO, w>eiwt dQ,<~>,



w(r, t) = ip B(t, 0) eiwt d W(r, w) =_% $, J,(r) GO, w) [email protected]).


As before we equate squared amplitudes, take expected values and assuming the existence of density functions we get the displacement evolutionary spectral density as .M, r, w) = 1% w)l fY(r, w) = g 2 f,(r)fs(r) GO, w) G,*(t, ~)_I&J),,. n=I s=l


Now using (79), (55) and (56) we have for the generalized force in the nth mode m


M,{G’,+ 2G,(i~ -+-5, CO”+) G,(w~ - w2+ 2i5, cm,,) eiw’dQ,b>>



f,(r,)e’“‘A(t,w)dP(r,,w)dr, =SI --m




As before we have the two representations for L,(t) as L,(t) = 1 Sf.(rJe’“’ --mA Choose such that


D,(t, w)/A(t, 0) = h,(w), and then restrict the random variables P(r, w) and Q,(w)

E{df’(r,, and

A(r, w) dP(r, w) dr, = f D,(r, CO)eiWtden(w).


w)> =fXr,,

r2, w)dw




The ensemble averages give the stationary spectra that we finally obtain, i.e. A(t,w) and G,(t, o) -+ unity. Since

G,,(t, w) -+ 1, t--tm

QU, w> -+ -J&



= zn

and so

h”(U) = Z”(W).

Then we obtain the differential equation G, + 2&iw + 5, WJ + G,(oz, - w2 + 2i5ow~) = A(t, w) (w;f - w2 + 2i5, ww,) from which we can obtain G,,. To determine the form offa(w thatfa(w).s =~~(cLJ)~~. Also _&(t,rl,r2,W)=


for the time-varying case we make use of our restriction



wheref,(t,r,,r,,w) is the evolutionary cross-spectral density of the excitation. So we can rewrite the expression forfz(w),s in (61) as

fn(rl)~(r2>f,(t,rl,r2)dr AI A2



Now substituting (83) into (80) we obtain the result for the evolutionary spectral density in terms of the excitation evolutionary spectral density as

f& r, ~1=

f.W_fXr)W, ~1G,*kw> IACt, w>12


’ 2.


Sn(r1)fs(r2)&(t,rlYr2,~)drI A1AZ

dr2. (84)

We must solve (82) to determinef,(t,r,w) just as for the single-degree-of-freedom For the particular case A(t, co) = 1 (stationary excitation), this reduces to (67).


7. AN EXAMPLE We have obtained equations (77) and (84) in the preceding sections describing the output evolutionary spectral density of a single- and multi-degree-of-freedom system, respectively, in response to an evolutionary input process. The equations relating the evolutionary spectra are identical to the standard equations relating stationary spectra except for the terms jG/A12 in (77) and G,G:/IAI 2 in (84). These terms are due to the excitation of the system by the time-varying character of the input. It would be useful to determine the magnitude of these terms in order to be able to assess the error likely to be incurred in omitting them from the equations altogether. It would then be possible, under certain conditions, to relate evolutionary spectra just as stationary spectra are related. The only difference in any response calculations would then be the estimation of the evolutionary input spectral density using the technique put forward by Priestley [4], to be used in the conventional response equations. It is to be noted that the results (77) and (84) were obtained by restricting the nonstationarity to a particular form. Also, they contain the excitation modulation function A(t, W)which is not calculated explicitly in the estimation of the evolutionary spectral density and consequently the equations would appear to be of academic interest only. However, in section 4 above we empirically derived an expression A(t,w) for a real process and in this section we shall go on to determine the effect of the term 1G/A I2 in a single-degree-of-freedom system for a particular situation. Despite the restrictions imposed on the form of nonstationarity the results are likely to give a general indication of the effect of the time-varying terms even though we will not be able to draw very detailed conclusions from the calculations.





For simplicity we restrict our example to a single-degree-of-freedom system. The modulating function is given by equation (27) and the differential equation to be solved is G + 2G(iw + SW,,)+ G(w; - 0J2+ 2i5WW~)= A(t, w) (0; - 0J2+ 2&Jwo). We recall that u(t) is the change in velocity and contains the time dependence. To solve this first substitute G = G, + iGt where G, and G1 are real. Equating real and imaginary parts of the resulting equation yields coupled second-order differential equations, and writing +=X2, :. x2=& G, = XI, and Gi = xj, :. 2, = xq, Gi = xq, gives the first order vector-matrix differential equation 0 2w 1 where Ai = A(wZ,- w2) and A, = A25ww0. The solution of this set of equations is obtained by using a Runge-Kutta technique which was coded in Fortran for the Southampton University I.C.T. 1909 Computer. 7.1.


Since the function A(t,w) was obtained by considering frequencies between 100 Hz and 1000 Hz this dictates the limits for the natural frequency we can choose for the second-order system. We shall choose an undamped natural frequency of 250 Hz. For the values of damping we shall consider, this will be very close to the resonant frequency of the system. This value of 250 Hz is also reasonably representative of an aircraft panel structural frequency, though perhaps on the high side. The values of damping we shall consider range from 5 = O$lO5to 5 = 0.1. Variations of IG/A12 are obtained for a frequency equal to the resonant frequency and also for frequencies some way below and above resonance, i.e. 200 and 300 Hz. The physical situation assumed is that the single-degree-of-freedom system is responding in a stationary manner to a stationary input for all t up to t = 0. Then for a time interval 0 < t G t, the excitation is time dependent and for t > t, is again representative of a stationary process identical to that represented by the record for t < 0. This means that the A(t,w) being used is unity for all t not between 0 and ta. Since the modulating function was obtained from jet noise considerations, this means that the engine is running with constant jet velocity tjO(=2000 ft/sec) for t < 0 and t > tp and varies in some prescribed manner for 0 < t < t,,. A convenient form of velocity variation will be assumed, namely, u(t) = a(1 - cospt), = 0



a is chosen as 150 ft/sec. Thus the effective jet velocity (I’,, - u(t)) has the appearance shown in Figure 4. The values of the period tD and of 5 will control the result so the effect of their variations are examined. Since the process is time-independent up to t = 0, the initial conditions for the differential equation are at t=O:

G, = 1.0 (=x,),

G, = 0 (=x,),

G;,= 0 (=x2),

G1 = 0 (=x4).



The solutions are obtained for the first half-period in all cases since the second half of the excitation time variation is merely a mathematical convenience.

Figure 4. Jet velocity profile. 7.2.



The results for 1G/AI* are shown graphically for a frequency equal to the resonant frequency in Figure 5 for different values of 5 for a time period of 10 set and in Figure 6 for varying t, with 5 = 0.005. Immediately apparent is the expected result that the lower the damping the greater is the deviation from unity of /G/Al”. But the deviation is very small even for 5 = 0.005 and for





tp= IO WC













Time (WC)

Figure 5. Effect of damping 5 on factor /G/Al*. I











Gi -




Time (WC)

Figure 6. Effect of decreasing tP on factor 1G/Al*.






t, = 10 set amounts to just over 4% at its largest. This is not surprising since even for 5 = 0.005, the natural frequency of 250 Hz means that the impulse response function of the system is short lived compared with the duration of the non-stationarity. The maximum value of 1G/A 12 is seen to occur just following the time when the rate of change of velocity is greatest (i.e. 4 period). Figure 6 is included to show the effect of increasing the severity of the time variation for a lightly damped system (5 = 0.005). We see that the error increases sharply as ta is reduced. We define a “characteristic time” of the system as the time taken for the envelope of the impulse function to reach l/20 its initial value. For a single-degree-of-freedom system this envelope is given by e-~~o* where time is measured from t = 0. It is normalized to equal unity at t = 0. For the system with 5 = 0.005 andf, = 250 Hz this characteristic time t, = O-381 sec. Taking t, as a characteristic time of the non-stationarity, then for a ratio t,/ts N 10 there is a maximum error of 11%. From Figure 6 we see that this error is increasing sharply as the ratio is reduced, so using conventional response equations for lightly damped systems, whose characteristic time t, is more than about a tenth of some reasonable characteristic time of the non-stationarity, is likely to be considerably in error. We must note here that since we have considered the first half period only the effect of the time-varying term is to augment the evolutionary spectral density; for the second half of the cycle the spectral function would be reduced by the same amount. However, we are concerned with “worst case” conditions and so will confine our attention to the first half cycle. The results for frequencies f = 200 and 300 Hz are not shown since the errors did not exceed 0.05 % even for 5 = 0.005. From these calculations we can conclude that the time-varying spectral density function f,(t, co) obtained from the response of a second order single-degree-of-freedom system whose resonant frequency is of the order of a typical aircraft panel structural frequency (and with similar damping) to a non-stationary excitation with spectrum &(t,w) is expressible by fy(t,ca) =fx(t,w)Ia(w)12 when the non-stationary trend is slowly varying compared with the “characteristic time” of the system (i.e., t,/ts > 10) for a maximum error whose absolute value is less than 11%. This discussion applies also to multimode systems. For a lightly damped system it has been shown [9] that the great portion of the output spectrum is contained in the direct terms (i.e., terms with n = s) and that the contribution to the total response of the cross terms (n # s) amounts to less than 10% for [ = 0.005. These calculations were carried out for a nine span beam with low damping in response to a convecting pressure field and can be considered as having the general characteristics of many practical multimode systems. Thus the response of a lightly damped multimodal system can be looked on as the sum of the responses in the individual modes, each of which is influenced by the non-stationarity as is the single-degree-of-freedom system. If, on the other hand, the cross acceptances are being included (which would be the case if modal damping were significantly higher), then the terms G, Gf/I A I2 (n # s) in (84) would need investigation though their effect is likely to be small owing to the increased damping. 8. CONCLUSIONS Spectral representations for non-stationary processes can, in certain cases, be very useful. The time dependent generalized spectra are extensions of the conventional spectral density functions. If the non-stationary process is known to be separable in the time domain as y(t) = m(t)Z(t) where Z(t) is stationary and m(t) contains the non-stationarity, the theoretical spectral function in section 2 above, the generalized spectral density function, can be used for investigations of the response of linear systems. 28



For a process known to be slowly time-varying, Priestley’s evolutionary spectral function can be constructed and in section 5 the non-stationarity response of a linear single-degree-offreedom system and multimode systems to a stationary excitation, using the evolutionary spectral density, is derived. The results agree identically with those of Caughey and Stumpf [5] and Lin 173who considered the same problems by other methods. Finally, by constructing a non-stationary excitation, the corresponding non-stationary response of a linear single-degree-of-freedom system is investigated by considering what effect the time dependent term due to system time-varying response has on the overall result. It is seen that even for a lightly damped, high frequency (5 = 0405, f = 250 Hz) system the error incurred in neglecting this term at the resonant frequency is less than about 11% for an “artificial”periodic non-stationarity with period four seconds. Thus it is reasonable to say that evolutionary spectra calculated using Priestley’s technique can be used directly in current response equations for typical air/spacecraft structures in a large number of practical situations where conditions change smoothly. For systems whose characteristic times are at all comparable with the non-stationary time trend (e.g. t,/ts -c 10 %) the error increases rapidly however, and in these instances time domain analysis is at present the only sure method of approach. ACKNOWLEDGMENTS

This research was carried out under a grant from the Science Research Council. The author is grateful to Professor B. L. Clarkson and Dr. C. A. Mercer for many helpful discussions in connection with this problem. REFERENCES 1. J. S. BENDAT and A. G. PIERSOL 1966 Measurement and Analysis of Random Data. New York:

John Wiley and Sons. 2. J. B. ROBERTS1965 J. Sound Vib. 2,3. On the harmonic analysis of evolutionary random vibration. 3. M. B. F%IESTLEY1967 J. Sound Vib. 6, 1. Power spectral analysis of non-stationary random processes. 4. M. B. PRIESTLEY1965 Jl. R. statist. Sot. B, 27, 204. Evolutionary spectra and non-stationary processes. 5. T. K. CAUGHEYand H. J. STUMPF1961 J. appl. Me&. Series E 28, No. 4. Transient response of a dynamic system under random excitation. 6. S. H. CIUNDALL @d.) 1963 Random Vibration. M.I.T. Press, Volumn 11, Chapter 3. 7. Y. K. LIN 1963 J. acoust. Sot. Am. 35,227. Non-stationary response of continuous structures to random loading. 8. A. POWELL1958 J. acoust. Sot. Am. 30, 1130. On the fatigue failure of structures due to vibrations excited by random pressure fields. 9. C. A. MERCER1965 J. Sound Vib. 2,3. Response of a multi-supported beam to a random pressure field.