IN. 1. NonLinear Mech‘uics, Vol. 25, No. 5. pp. 493509. 1990 Printed irl Great Britain.
CO2&7462,W S3.00 + .W Pergaman Press pk
RESPONSE STATISTICS OF NONLINEAR SYSTEMS TO COMBINED DETERMINISTIC AND RANDOM EXCITATIONS A. H.
NAYFEH
and S. J.
SERHAN
Engineering Science and Mechanics Department, Virginia Polytechnic Institute and State University, Blacksburg, VA 24061, U.S.A. (Received 17 June 1988) AbstractA secondorder closure method is presented for determining the response of nonlinear systems to random excitations. The excitation is taken to be the sum of a deterministic harmonic component and a random component. The latter may be white noise or harmonic with separable nonstationary random amplitude and phase. The method of multiple scales is used to determine the equations describing the modulation of the amplitude and phase. Neglecting the thirdorder central moments, we use these equations to determine the stationary mean and meansquare response. The effect of the system parameters on the response statistics is investigated. The presence of the nonlinearity causes multivalued regions where more than one meansquare value of the response is possible. The local stability of the stationary mean and meansquare responses is analysed. Alternatively, assuming the random component of the response to be small compared with the mean response, we determine steadystate periodic responses to the deterministic part of the excitation. The effect of the random part of the excitation on the stable periodic responses is analysed as a perturbation and a closedform expression for the meansquare response is obtained. Away from the transition zone separating stable and unstable periodic responses, the results of these two ap proaches are in good agreement. Comparisons of the results of these methods with that obtained by the method of equivalent linearization are presented.
1. INTRODUCTION
Most of physical excitations exhibit randomly fluctuating characters and contain a wide spectrum of frequencies that may result in severe vibrations. Due to the highly unpredictable trends of these natural hazard excitations, a stochastic formulation is needed to describe them. They can be completely characterized by probability density functions or certain statistical measures. This explains the wide range of applications of random vibrations to the solution of practical engineering problems. The general problem of random vibrations has been investigated in the literature for many years through the study of Brownian motion. Einstein [l] was the first to study the problem, which was later extended and developed by Fokker [2], Planck [3], Kolmogorov [4], Feller [YJ, Ito [6] and Doob [7]. Over the years, a number of different techniques has been used for the study of nonlinear mechanical and structural systems disturbed by random excitations. So far, exact solutions of the FokkerPlanckKolmogorov equations are only available for a very limited number of problems. Thus, approximate methods have been developed and used to treat many of these problems. These include the method of equivalent or stochastic linearization, perturbation methods, stochastic averaging and series expansions. Dynamical systems subject to random excitations exhibit nonlinear behavior for sufficiently large motion [8,9]. The sources of the nonlinearities may be material (nonlinear elasticity, plasticity, strain hardening, etc.), geometric (excessive displacements, curvature, etc.), damping (dry friction, hysteretic damping, etc.), and inertia. When a linear system is excited by a Gaussian input, the output will be Gaussian. However, the output of a nonlinear system to a Gaussian input is not Gaussian except in some special cases [lo]. The deviated output yields different probability functions that influence the investigation of the failure criteria of systems [11131. Linear problems can be solved exactly by a variety of techniques [14,15]. A comprehensive treatment of linear systems subject to random excitations is available in standard textbooks such as those of Clough and Penzien [16], Crandall and Mark [17], and Lin [18]. For a comprehensive survey of the recent Contributed by P. D. Spanos. 493
A. H. NAYFEH and S. J. SERHAN
494
developments in the area of nonlinear random vibrations, the reader is referred to Caughey [19]. Ibrahim and Roberts [20], Spanos [21], Roberts [22241, Crandall and Zhu [25], To [26], and Roberts and Spanos [27]. In this paper we describe a secondorder closure method based on the method of multiple scales to determine the mean and meansquare response of nonlinear singledegreeoffreedom systems to an excitation that is the sum of a deterministic primary resonant component and a random component. The method is described using the DuffingRayleigh oscillator ii + c$u  ES1ti + +.s6,ri3 + &au3= E[hOcosRt  k, sin!& + t(t)]. (1.1) Here, < is a general zeromean random excitation. The excitation is taken to be O(E)so that the primary resonance, the damping and the nonlinearity balance each other. In the present analysis, the random component of the excitation is taken to be either white noise or sinusoid having a random amplitude and phase. The method of multiple scales [28,29] is used to obtain two averaged firstorder differential equations. They describe the modulation of the amplitude and phase of the response with the slow time scale. Three alternate approaches are used to determine the mean and meansquare values of the amplitude and phase from these modulation equations.
2. THE METHOD OF MULTIPLE SCALES To use the method of multiple scales [28,29] for analysing the case of primary resonance (i.e. Q = oe), we describe the nearness of the resonance by introducing the detuning parameter v defined according to R = wg + EV. (2.1) Then, we seek a uniform approximate
solution of equation (1.1) in the form
I.4= u,(T,, 2i) + EUi(&, r,) + . . .
(2.2)
where To = t is a fast scale, characterizing motions occurring with the frequencies o,, and R, and T, = Et is a slow scale, characterizing the modulation of the amplitude and phase with the nonlinearity, damping and resonance. In terms of the T,, we use the chain rule to transform the ordinarytime derivatives into partial derivatives according to d =DD,+&D,+ dt c
dt2
...
(2.3)
= D2 + 2&D01D + ... ’
where Dn = d/aT,. Substituting equations (2.2H2.4) into equation coefficients of like powers of E to zero gives Order co:
D;u,
(1.1) and equating
+ co;uo = 0.
(2.5)
Order E: D;u,
+ co&,
= 2D,D,u,
 f62(Douo)3
+ 6,D,u,
 au; + h o cos QT,  k, sin LIT, + 5.
(2.6)
The general solution of equation (2.5) can be written as u. = A(T, )eiooTo + cc
(2.7)
where cc stands for the complex conjugate of the preceding terms. To this order, the function A is arbitrary and it is determined by imposing the solvability condition at the next level of approximation. Substituting equation (2.7) into equation (2.6) yields Diu, + co~ul =  io,(2A’  6,A + ~~~2A2A)eiooTo
+
+jb2,,;A3e3iooTo
+ $,e'"T"
+
_
jikoeiRTo
a(A3e3iwoTo
+ cc + 5.
+
3A2AeimTo)
(2.8)
495
Response statistics of nonlinear systems
Using equation (2.1) to eliminate the terms that lead to secular terms from equation (2.8) yields + 3aA2A  $hoeiVTI  $ikOeiYTL
io,(ZA’  6,A + oi6,A2Z)
(2.9)
Starting from the averaged equation (2.9), we used three methods to determine the statistics of the response. 3. LINEARIZATION
METHOD
A lowintensity noise is assumed to separate the strong mean motion from its weak fluctuations. A similar method was used by Berstein [30] and Stratonovich [31] to determine the influence of noise on the amplitude of the limit cycle of the van der Pol oscillator and by Nayfeh [32] to investigate the response of the van der Pol oscillator with delayed amplitude limiting to a sinusoidal excitation. Nayfeh chose the amplitude of the excitation to be either constant or white noise. For the latter case, he obtained the conditional probability function for deviations of the amplitude from its limit cycle value. It is convenient to express A in the polar form A = &rei(vrI Y)_
(3.1)
Substituting equation (3.1) into equation (2.9), separating real and imaginary parts, and replacing Ti with EC,we obtain ci =&2(6,
*w$6,a’)+$siny
+$coSy
+ 5i
(3.2)
c2 = 2n tcos4di$, 2nwo s 0
(3.4a)
0
3&a aj,=~vaa3+oCOSy
8~0,
.sh
0
Eke
Ksiny+r2
2~0,
0
where E
f#J= Rt  y.
(3.4b)
According to this approach, we assume that the C;,,are small compared with ho and k,. Next, we determine the steadystate values a, and y. of the amplitude and phase when
ko
0
3a &a
ho :  vao =zw,cosyo
(3Sa)
0
ko my,.
.
(3Sb)
Squaring and adding equations (3.5) yields the frequencyresponse w6b2 0 2 + 9a2 160; a8 
w$5,6,
+
20,
6av
k2
equation h;
u; + (6: + 4v2)&  $  7
= 0.
(3.6)
WO
Hence, to the first approximation,
the deterministic steadystate response is given by
u=a,cos(Rty,)+
...
(3.7)
where a, and y. are given by equations (3.5) and (3.6). Next, we determine the effect of the noise (i.e. t1 and r2 # 0) on the deterministic steadystate motion. To this end, we let a=a,+a,, Substituting
y=yo+y1.
(3.8)
equations (3.8) into equations (3.2) and (3.3) and linearizing the resulting
A. H. NAYFEH and S. J. SERHAN
496
equations, we obtain
Cx(t) + G(r)
k(t) =
(3.9)
where (3.10)
G(t)={;;;:;}. C=[;;]
(3.11)
520,Yo) s1w = 51k Yo), g*(t) = ~
(3.12)
a0
EV
cj =  u a0
9&a
o
and
8wo
cq = .)s(6i  &&,aij).
(3.13)
We note that C is a real nonsymmetric matrix. The general solution of equation (3.9) can be obtained by modal analysis. The eigenvalues and eigenvectors of C are given by cx = Lx
(3.14)
while the eigenvectors of the adjoint problem are given by CTV = Iv.
(3.15)
Then, the modal matrices are
a,=cIc3[czc3+(1,cr)*]i
fern=
It is convenient to introduce a transformation
(3.17)
1 and2.
that decouples equations (3.9). To this end,
we let
(3.18)
x(c) = X2(t). Substituting equation (3.18) into equation (3.9) and premultiplying
by VT yields
21(r) = JiZl(Q +
91
(3.19)
i*(t)
v2
(3.20)
=
~*z*(Q
+
where tf,=r1+
1 et,
for n = 1 and 2.
Then, the statistics of the uncoupled coordinates
z1 and z2 are
(zwl> = 0 (Z,(r1)Z,(t2))
= e’“&+“““’
(3.21)
(3.22a)
I’ r* e&r1 +ZnQ). ss0 0
form, n = 1, 2
where the angle brackets stand for the expectation. As examples, we consider the cases of sinusoids having nonstationary tudes and white noise.
(3.22b) random ampli
Response statistics of nor&near
3.1, constatianar~randomamplitude
systems
497
sinusoids
In this case, we Iet if(t) = h, cosRt  k, sinRt
where h, and k, are separable nonstationary hi =f,(T,)w,,
(3.23)
random functions described by (3.24)
k, =.Mq)wz.
Here, thef, are slowly varying deterministic functions of time, whereas the w, are slowly varying stationary random functions. Road roughness and seismic ground motion could be modeled by separable nonstationary random processes. Performing the integrations in equations (3.4a) and using equations (3.21), we can write the q,, as q.(t) = d,,h,(t)
 d,,k,(t)
4l Cl cm Yo aOc3 1 1
where d
4, = &
for n = 1, 2
A, 

Cl
cosy, t siny,
.
(3.25)
(3.26a)
(3.26b)
aoc3
0
A narrowband random excitation is a special case of the excitation given by equations (3.23) and (3.24). Here, t(t) is chosen to be a zeromean Gaussian narrowband random excitation. It could be obtained by filtering a white noise through a linear filter [31,33,34j; that is, z t n25 + 84 = /?‘%2 W (3.27) where !J is the center frequency of 5, /l is the bandwidth of the filter, the autocorrelation function of the white noise W is given by (3.28)
R,(r) = 2?rS,6(r)
and 6 is the Dirac delta function. Substituting equation (3.23) into equation (3.27) and performing deterministic and stochastic averaging of the equations describing the modulations of hi and k, with time, one obtains (3.29a)
rt,[email protected],
= [email protected]
W,.
(3.29b)
The white noise components W, and W2 are independent and their autocorrelation functions are given by equation (3.28). The autocorrelation functions of h, and k, are R,,,(Z) = R,,(T) = 7cSoe81’112.
(3.30)
The correlation time of h, and k, is 0 (l/j?). This means that for sufficiently small values of 8, h, and k, are slowly varying functions of time. Using equations (3.22b), (3.25), (3.26) and (3.30), we obtain the following stationary statistics of the uncoupled coordinates z,:
To the first approximation,
the solution of equation (1.1) can be written as
u = a, cos (Qt  yo) + q cos(Qt  yo) + yt a, sin (S2t  yo) + . . .
(3.32)
where o. is given by equation (3.6), y. is given by equation (3,5a), and the statistics of a1 and y1 are &>=O, (a:>
= a:+:>
+ ai
=O + 2a,a,(z,z,)
(3.33a) (3.33b) (3*33,_)
498
A. H. NAVFEH and S. J. SERHAN
(3.33d) Using equations (3.8) and (3.33), we obtain the following statistics for the amplitude and phase of the response: (a) = a,, is; = Of, (3.34) (‘i> = 70,
(3.35)
a,” = &
3.2. White noise In this case, we have R&r) = 27&6(r).
(3.36)
The autocorrelation function given by equation (3.36) implies that < is a fast varying random function of time. Terms including t(t) in equations (3.2) and (3.3) can be replaced by their means and fluctuating components [31]. This stochastic averaging is restricted to random excitations with small correlation time compared with the relaxation time of the oscillator. The main advantage of stochastic averaging is to decouple the amplitude and phase equations (3.2) and (3.3). However, the presence of the dete~inistic excitation will keep this coupling. It is applicable to any lowintensity random excitation. Using equations (3.4), (3.21), (3.22), and (3.36), we can express the expected values and the meansquare values of the uncoupled coordinates z, as (3.37a)
= 0
cw7l> = 
e= ns, w$.,
1+(’ 24
+ i.,)
SECONDORDER
4.

clfazl Cl) a$$
CLOSURE
for n, m = 1,2.
(3.37b)
1
METHOD
The dete~inistic motion given by equations (3.5) and (3.6) has no feedback from the perturbations due to the noise. This limits the applicability of the linearization method in the transition zone separating stable and unstable periodic responses. In fact, the preceding linearization scheme predicts infinite meansquare response at the bifurcation points separating stable from unstable deterministic responses. To overcome this problem, we propose a secondorder closure method. To describe the method, we find it convenient to express A in equation (2.9) in the form A = 4(x  iZ)eivTI.
(4.1)
Substituting equation (4.1) into equation (2.9), separating real and imaginary parts, and replacing Ti with Et yields k=

+3,X
i = jES,Z
eke 
E
?jsin Rt dt
(4.2a)
Znloo < cos Rt dt + &V,X+ &ho + 5 2wo 2lr s ()
(4.2b)
svez
+ 2wo
28
where (4.3a)
6c =  s, t *0&5,(x2 + 22) v, = v 
$(x2+22).
(4.3b)
0
Next, we separate each of x and z into mean components components y, and y, according to x = xg + y,
and
z=z,+yZ.
x0 and z. and random (4.4)
Substituting equations (4.4) into equations (4.2) and (4.3) and separating mean and fluctuating components, we have:
Response statistics of nonlinear systems
499
mean motion 1 1 1 20 = j&6,X0  [email protected]*x;  &o;b,x,z; 8 8 3&a  GZO 0
3E g~;wo
+ +
[
3.sa  EVZO + x;zo 8~0
1 [
1
 p&%X0
3&a + =z;
1
0
(4.5a)
j$?
0
1 1 zo = &SIZo  &o;b,x;z, 2 8
1 3&a  [email protected]*z; + &VXO  x; 8 80,
36a  xoz; 8wo
(4.5b) fluctuating motion f(t) = Cy(t) + G(t) + WY,, yz),
(4.6)
Znloo 910)
=

&
(4.7)
s 0
1 3 Cl = T&6,  j&CQ2x;
1  @6,Z;
3&a I xozo
(4.8a)
4WO c2 =

c3 =

1 &w;b,x,z,
4
1
&o;b,x,z,
4
3Ea
EV+GX +
EV 
9&a
x;
80,
9Ea
; +Gz;
(4.8b)
3&a z;
(4.8~)

SW,
1 1 3 3&a cq = T&6,  8&W;6ZX;  B&o;6,z;  xozo
(4.8d)
4WO
?ll
=
(Y:>)+xo(Y; (Yf>)

+ 2Zo(Y,Y,  (Y,Y,))+Y: +YlY:  (YlY:>l +~C3zo(Y:  )+ 2XOcylY2  (Y1Y2)) 0
Y: (vi>+Y:Y2  (Y:Y,)l 1 n2= p&%CzO(Y:  )+2xocylY2 
(4.9a)
+ 3Zo(Yi  (Y:N+Y:Y2  (Y:Yd +y: (Y:N gc3xo(Y:  )+2zo(Y1Y2 
Y: +Y,Yt 
(4.9b)
A. H. NAYFEH and S. J. SERHAH
500
For noise whose intensity is small compared with h, and k,, y, and y2 are small compared with x0 and z,,, and hence central moments higher than the second can be neglected in equations (4.5) and (4.6). Consequently, it follows from equations (4.5) that the steadystate mean motion is given by 1 1 1 &Six0  EC&&x;  s0+,x,z; 2 8 8
3&Z + z; 8~0
 &VZO 3&X + x&) 8%
(4.10a) 1 2
E81Z0

1 8
EW;d,X;Z,

1 8
&CO&Z;
+
3&a EVXO 
Xi
3E2 
80,
XoZ;
80,

(4. lob)
We neglect central moments beyond the second, and hence neglect N in equations (4.6). Using modal analysis to decouple these equations, we find that the statistics of the uncoupled coordinates zi and z2 for a narrowband random excitation are given by equations (3.22a) and (3.31) where d,,
=
Et4  Cl)
2woc3
, 4, = 
&0
(4.11)
For the case of a white noise, these statistics are given by equations (3.37). The meansquare values of y, and y, have the form described in equations (3.33) for a, and yl. To the first approximation, the solution of equation (1.1) can be expressed as u =(x0 + y,)cosRt
+ (z. + y,)sinRt
+ ...
(4.12)
Equations (3.33) and equations (4.10) and either equations (3.31) or equations (3.37b) represent a system of five simultaneous nonlinear algebraic equations. A modified Powell’s hybrid algorithm is used to solve this system. For the case of a narrowband random excitation, premultiplying equations (4.6) and (3.29) by y,, y,, h, and k, and taking the expectations of the resulting equations, we have (4.13)
P=BP+Q where P = {, (~9,
(Y~Y,),
(y,h,)v
(~,k,)j~
(4.14)
and B and Q are given in Appendix A. The stability of a particular fixed point of equations (4.10) and (4.13) to a perturbation proportional to exp(It) depends on the eigenvalues of C and B. 5.
SEMILINEARIZATION
METHOD
Rajan and Davies [34] developed a method for the investigation of the response of the Duffing oscillator to a zeromean narrowband random excitation. In this section, we generalize this method for the case of nonzeromean random excitations and point out its limitations. To determine the stationary meansquare value (x2 + z2) of the amplitude, we follow Rajan and Davies and assume average values for 6, and v,. Substituting equation (3.23) into the integrands in equations (4.2) and performing the deterministic averaging, we find that the integral terms in equations (4.2a) and (4.2b) become Ekl /2w, and ch1/2w0, respectively.
Response statistics of nonlinear systems
501
Then, we combine equations (4.2) into
d(xf + z2>
= b&x2
dt
+ t=) + $(zh
(5l)
I xk)
where h = ho + h,
and
k = k, f ki.
(5.2)
Next, we add k times equation (4.2a) to h times equation (4.2b), use equations (3.29), and obtain d(zh + xk) = v,(xh  zk)  ;(& + jT)(rh + xk) dt f
;Pilt,(z) +koCx>)+ {h2+kz>+
gzw,
2o
+ XW,).
J
0
(5.3)
Subtracting k times equation (4.2b) from h times equation (4.2a) and using equations (3.29), we obtain d(xh
 zk) dt
=
 v, (zh t xk)  ;(& + Bl(xh
+
 zk> f fS(ho(x>
J;
zw,>.
 ko(z>)
(5.4)
Taking the expected values of equations (4.2) yields
d = dt
W> = dc

+x>
20
(5.5)

fd,(z>+&V,(X) + 2.0
(5.6)
Eve(Z) +
Following Rajan and Davies, we assume that (x Wi) = 0 and (zW,) = 0. Then, the stationary solutions (long time averaging) of equations (5.1) and (5.3H5.6) are given by
where Fi = +(hs + ki)
and
( F2) = Fi t +(hf + kf).
(5.8)
A shorttime averaging is implied in equation (5.8). Equation (5.7) is a nonlinear algebraic equation in (x2 + z2). The response is assumed to be pseudo sinusoidal and the frequency shift and damping parameters of the system are given by s c = 8,
+ .&f#,(x*
4 22)
(5.9)
v,=v~(x2+z2). 0 Richard and Anand [35] used a generalized form of van der Pal’s method to investigate the stability of the response of the Duffing osciitator to a narrowband random excitation. With probability one, they found the conditions under which the perturbations of the averaged equations tend to zero. Davies and Nandlall [36] used the phaseplane diagram [(r.i* ) vs (u* )] for their stability analysis. Meansquare values are obtained from the time dependent FokkerPlanckKolmogorov equation. A Gaussian closure approximation is adopted to handle higherorder moments. Fast varying terms are averaged out to obtain smooth time histories of (uz> and (ti*). The stability of a particular fixed point of equations (5.1). (5.3)(5.6), (5.9) and (5.10) to a perturbation proportional to exp(lt) depends on the eigenvalues of the Jacobian matrix.
A. H.
502
NAYFEH and S. 6.
J.
SERHAN
RESULTS
In the pvplane, where p = l4w0u0, 2 z the frequencyresponse equation (3.6) provides a oneparameter family of response curves. When ‘1= 0, the response curves are symmetric with respect to the paxis. Let fo = (hi t k,) * li2 . Whenf, = 0, the curves of the family degenerate into the line p = 0 and the point v = 0 and p = 1. As fo increases, the curves first consist of two branchesa branch running near the vaxis and a branch consisting of an oval, which can be approximated by an ellipse having its center at (0, 1). As f0 increases, the ovals expand and the branch near the vaxis moves away from this axis. When f$ reaches the critical value 16/27 a2(6 = 6, = a,), the two branches coalesce, and the resultant curve has a double point at (0, j), as shown in Fig. 1. As f0 increases further, the response curves become open curves, which continue to be multivalued functions until fi exceeds the second critical value 32/27 6’. Beyond this critical value, the response curves are singlevalued functions for all values of v. When a # 0, the response curves lose their symmetry with respect to the p axis and bend to the right for a > 0 (hardening nonlinearity) and bend to the left for a < 0 (softening nonlinearity). Representative curves are shown in Fig. 2 when a > 0. When fo = 0, the family of
1 .co Fig. 1. F~quencyresponse
0.67
0.33
0.00
0.33
0.67
.oo
curves For primary resonances of the Rayleigh oscillator for 8 = 1.0, w0 = 1.0, k, = 0.0 and h, = k, = 0.0.
UNSTA3LE
1.02
3.67
0.33
0.00
0.33
0.67
I. 00
Fig. 2. Frequencyresponse curves for primary resonances of the DuffingRayleigh oscillator for (r = 1.0. z = 0.2. co0 = 1.0. k, = 0.0 and h, = I;, = 0.0.
Response statistics of nonlinear systems
503
curves degenerates into the line p = 0 and the point p = 1 and v = 301/20& As f0 increases, the response curves consist of two branchesa branch running near the v axis and a branch consisting of an oval. As f0 increases further, the oval expands and the branch near the vaxis moves away from it. When fo reaches a critical value fat, the two branches coalesce and have a double point. As fo increases beyond fat, the corresponding response curve continues to be a multivalued function of v until a second critical value fo2 is exceeded. For all values of f0 greater than &, the response curve is a singlevalued function of v. Not all of the solutions given by the frequencyresponse equation (3.6) are realizable because some of them are unstable. The solution given by the frequencyresponse equation (3.6) is stable if none of the eigenvalues defined in equation (3.14) has a positive real part. The eigenvalues are given by the characteristic equation ,I2  e&l  2p)L + &‘A = 0
(6.1)
where 6vap
A = $S2(1  4p + 3p2) + v2  
+ 27a2p2
___ 4w;
4
and
1 p = [email protected];.
(6.2)
When A c 0, the roots of equation (6.1) are real and have different signs. Hence, the predicted deterministic steadystate motions correspond to saddle points and are unstable. These correspond to the interior points of the ellipse A = 0 in the pv plane. Since the discriminant of equation (6.1) is 24vup 27u2p2 D = g2p2  4v2 + 3 6
(6.3)
WO
00
periodic motions corresponding to D c 0 are focal points, while those corresponding to D > 0 are nodal points. For 6 > 0, these points are stable or unstable depending on whether p is greater or less than +. The dark lines in Figs 1 and 2 separate stable from unstable solutions; all solutions corresponding to points above the dark lines are stable and those corresponding to points below the dark lines are unstable. Figure 3 shows that the deterministic response predicted by the perturbation expansion is in good agreement with that obtained by numerically integrating equation (1.1) using a fifthorder RungeKuttaVerner technique (numerical simulation). The effect of the excitation frequency on the primary response is shown in Fig. 4. The linearization method is used in Fig. 4(a) to predict the effect of a narrowband random excitation on the deterministic response. Near the saddlenode bifurcation point (the fixed points lose stability with a real eigenvalue crossing the imaginary axis) and Hopf bifurcation points (the fixed points lose their stability with the real part of a complex conjugate pair of eigenvalues changing sign from negative to positive), the deterministic motion is sensitive to any outside perturbations and the linearization method fails. In Fig. 4(b), the secondorder closure method is used to obtain the meansquare response. The secondorder closure method provides the necessary feedback from the pertubations caused by the noise to the mean motion. This limits the secondorder moments predicted by the linearization method at the bifurcation points separating stable from unstable responses. Since a > 0, the response curves in Fig. 4 are bent to the right, causing multivalued regions. The multivaluedness is responsible for a jump phenomenon at the saddlenode bifurcation points that leads to an abrupt change in the response. The semilinearization method is based on assuming an average value (x2 + z2) of the amplitude. With this assumption, substituting equations (4.4) into equations (4.2) and (4.3) and separating mean and fluctuating components for the case of a narrowband random excitation, we have: mean motion i.
1 = j&6,X,
1 1  &w;b2X~  ~&w~G,x,z~ 
8 
1 3Ea &W$52X0 z8w 8
I 1 [
0 0
EVZO
g~W~~,Xo
+
3&a 3&a x;zo + z;
800
800
3ea 
GZO 0
1 (YS>
&ko
+ jg
0
A. H.
504
NAYFEH and
S. J.
SERHAN (a)
P) I
400
417
433
450
467
483
W
i 417
3co
433
450
467
483
500
t
Fig. 3. Comparison of the mean solution predicted by the perturbation analysis with that obtained by numerical simulation for 6 = 1.0. IX= 0.2, o,, = 1.0, h, = 1.414, E = 0.08, k, = 0.0, and h, = k, = 0.0: (a) numerical solution; (b) perturbation solution.
LINERRIZRTION A N cu, a
0 0
0.10
0.06
0.02
0.02
0.06
Y
2NDJJROER CLOSURE A N su, d
Fig. 4. Frequencyresponse curves (narrow band) for 6, = 0.02, 6, = 0.08. z = 0.1, 00 = 1.0, k, = 0.0316, h, = 0.0316, nS, = 0.00005, and /3 = 0.01: (a) linearization method; (b) secondorder closure method.
505
Response statistics of nonlinear systems
1
i,
2
= EblZg

1
3Ea
1
[email protected],X;Z, 8

8
&W&Z;
+
&VXO 
3ECl
80,  Gxoz’
Xi
(6.4b) fluctuating motion
(6.5a) 3Ea
3Ea
x ; z 8wo 80,
1 1 [E& 2
; Y, +
 ;EW;bZZ;
&0&X;
1
J’, + 2. 0
(6Sb) The mean motion given by equations (6.4) provides a feedback from the noise through ( y: ) and (y:). However, the cross correlation (yr y, ) has no contribution and the coefficients of the terms x0 (y:) and z,(y:) are onethird of those given by equations (4.5) obtained by the secondorder closure method. The linearization method and the secondorder closure method have the same form with respect to the fluctuating motion. According to equations (4.6), the fluctuating motion given by equations (6.5) does not include the xozo terms. The coefficients of xiy, and ziy, are onethird those given by equations (4.6). As a result, the accuracy of the semilinearization method is somewhere between the linearization method and the secondorder closure method. Figure 5 shows a comparison among the results of these three methods and that obtained by the method of equivalent linearization (details in Appendix B). Away from the transition zone between stable and unstable responses, the results are in good agreement. The linearization method is simple and direct but it fails to give good approximations to the meansquare response at the stability boundaries. The secondorder closure method is expected to be the most accurate but it requires the solution of a system of five simultaneous nonlinear algebraic equations. The local stability of the stationary mean and meansquare response is analysed. Figure 6 shows the effect of a low intensity noise, around the natural frequency o. of the system, on the stability of the meansquare response. The strength of the noise is taken as l/450 of the strength of the deterministic excitation. Since a = 0, the response curve is symmetric with respect to the (u’)axis. The response curve has its peak at the perfectlytuned frequency 0 = w,(i.e. v = 0). The mean motion loses stability when [VI2 0.093. The presence of the noise expands the stability region Iv12 0.1. As the amplitude of the mean excitation increases, the effect of the noise on the steadystate mean response decreases as shown in Fig. 7.
+
m 04 d
LINERRIZRTION
+
ZNDORDER

SEMILINERRIZRTION EO. LIN.
CLOSURE
I GRUSSIRNI
: d 0.10
0.06
0.02 Y
0.02
0.06
I. IO
Fig. 5. Comparison of the variations of the meansquare responses with the narrowband excitation frequency obtained using the different methods: 6, = 0.2, s, = o.z, a = 0.0, o0 = 1.0, h, = 0.2, k, = 0.001, nS, = 0.002 and B = 0.01.
A.
506
H.
NAYFEH and
S. J. SERHAN (a)
I
1
0
0.
IO
0.06
0.02 Y
0.02
0.06
0.10
STABLE L
I (‘4
0.
IO
0.06
0.02 Y STRBLE
L
0.02
0.06
0. IO
I
Fig. 6. Frequencyresponse curves for 6, = 0.2, 6, = 0.2, a = 0.0, o0 = 1.0, h, = 0.3, k, = 0.001, nS, = 0.0001 and B = 0.01: (a) no noise; (b) with noise.
1
01
0.03
0.09
0.16
0.22
0.29
0.35
0.42
ho Fig. 7. Effect of the mean amplitude of excitation upon the meansquare response to white noise for 6, = 0.02, 6, = 0.08, 01= 0.05, o,, = 1.0, k, = 0.0, nS, = 0.005 and v = 0.03.
7. CONCLUSIONS
A secondorder closure method is presented for determining the response of nonlinear systems to random excitations. As an example, the method is used to determine the response of the DuffingRayleigh oscillator to an excitation consisting of the sum of a deterministic harmonic component and a random component. The random component is assumed to be white noise or sinusoidal whose amplitude and phase have two partsa constant mean and a separable nonstationary random process. The method of multiple scales is used to determine the modulation of the amplitude and phase with the nonlinearity and the excitation. Three methods have been used to determine the mean and meansquare response from the modulation equations. The local stability of the stationary mean and meansquare response is analysed. For some range of parameters, the presence of noise can expand the stability region of the mean response. Comparisons of the results of these methods with the results obtained by the method of equivalent linearization are presented. The present analysis attacks the problem in a way that the statistics of the response can be determined for a general type of excitation and nonlinearity. The nonlinearity bends the frequencyresponse curves, causing multivalued regions in the meansquare response. The multivaluedness is responsible for a jump phenomenon.
Response statistics of nonlinear systems
507
AcknowledgementsThis work was supported by the Air Force Office of Scientific Research under Grant # AFOSR860090 and the National Science Foundation under Grant # MSM8521748.
REFERENCES 1. A. Einstein, Uber dievan der molekularkineteschen theorii der warme gefordert bewegung von in ruhenden flflssigkeiten. Ann. Phys. (Leipzig) 17, 549 (1905). 2. A, P.Fokker, Die mittlere energie rotierender electricher dipole im stahlungsfeld. Annals. Phys. 43,810 (1914). 3. M. Planck. Uber einen Satz der Statischen Dvnamik und Seine Erweiterung in der QuantenTheorie. Stizungsber. Preuss. Akadamie Weiss, 324 (1971). . 4. A. Kolmogorov, Uber die analytischen methoden in wahrsheinlichkeitsrechnung. Mathematische Annalen 1’04, 415 (1931). 5. W. Feller, Zur theorie der stochastischen
6. 7. 8. 9. 10. 11. 12. 13. 14. 15. 16. 17. 18. 19. 20. 21. 22.
prozesse (existenz und eindeutigkeitssatze). Marh. Ann. 113, 133 (1936). K. Ito, On stochastic differential equations. Mem. Am. Math. Sot., No. 4 (1951). J. L. Doob, Stochastic Processes. Wiley, New York (1953). R. H. Lyon, Empirical evidence for nonlinearity and directions for future work. J. Acousr. Sot. Am, 35, 1712 (1963). S. H. Crandall, Nonlinearities in structural dynamics. Shock Vib. Dig. 6, 2 (1974). Y. Yong and Y. K. Lin, Exact stationaryresponse solution for second order nonlinear systems under parametric and external Whitenoise excitations, J. appt. Mech. 54, 414 (1987). J. B. Roberts, First passage probability for nonlinear oscillators. ASCE J. Engng Mech. Div. 102, 851 (1976). J. B. Roberts, First passage time for oscillators with nonlinear damping. J. appl. Mech. 45, 175 (1978). J. B. Roberts, First passage time for oscillators with nonlinear restoring forces. J. Sound Vib. 56, 71 (1978). M. C. Wang and G. E. Uhlenbeck, On the theory of Brownian motion II. Rev. Mod. Phys. 17, 323 (1945). A. T. BharuchaReid, Elements of the Theoryof Markov Processes and Their Applications. McGrawHill, New York (1960). R. D. Clough and J. Penzien, Dynamics of Structures. McGrawHill, New York (1975). S. H. Crandall and W. D. Mark, Random Vibration in Mechanical Systems. Academic Press, New York (1963). Y. K. Lin, Probabilistic Theory of Structura! Dynamics. McGrawHill, New York (1973). T. K. Caughey, Nonlinear theory of random vibrations. Ado. appt. Mech. 11, 209 (1971). R. A. Ibrahim and J. W. Roberts, Parametric vibration; part V: stochastic problems. Shock Vib. Dig. 10, 17 (1978). P.T. D. Spanos, Stochastic linearization in structural dynamics. Appl. Mech. Rev. 34, 1 (1981). J. B. Roberts, Response of nonlinear mechanical systems to random excitation; part I: Markov methods.
Shock Vib. Dig. 13, 17 (1981). 23. J. B. Roberts, Response of nonlinear mechanical systems to random excitation; part II: equivalent linearization and other methods. Shock Vib. Dig. 13, 15 (1981). 24. J. B. Roberts, Techniques for nonlinear random vibration problems. Shock Vib. Dig. 16, 3 (1984). 25. S. H. Crandall and W. Q. Zhu, Random vibration: a survey of recent developments. J. appl. Mech. 50, 953
(1983). 26. C. W. S. To, The response of nonlinear structures to random excitation. Shock Vib. Dig. 16, 13 (1984). 27. J. B. Roberts and P.T. D. Spanos, Stochastic averaging: an approximate method of solving random vibration problems. Inr. J. NonLinenr Mech. 21, 111 (1986). 28. A. H. Nayfeh, Perturbation Methods. Wiley, New York (1973). 29. A. H. Nayfeh, Introduction to Perturbation Techniques. Wiley, New York (1981). 30. I. L. Berstein, On fluctuations in the neighborhood of periodic motion of an autooscillating system. C. R. Acad. Sci. U.S.S.R. Mot.Fiz 20 (1938). 31. R. L. Stratonovich, Topics in the Theory of Random Noise. Gordon and Breach, New York (1983). 32. A. H. Nayfeh, Forced oscillations of the van der Pol oscillator with delayed amplitude limiting. IEEE Trans. Circuit Theory CT15, 192 (1968). 33. S. Rajan, Random superharmonic and subharmonic response of a Duffing oscillator. Ph.D. Thesis, University of New Brunswick, New Brunswick, Canada (1987). 34. S. Rajan and H. G. Davies, Multiple time scaling of the response of a Duffing oscillator to narrow band random excitation. J. Sound Vib. (1987). 35. K. Richard and G. V. Anand, Nonlinear resonance in strings under narrow band random excitation, Part 1: planar response and stability. J. Sound Vib. 86, 85 (1983). 36. H. G. Davies and D. Nandlall, Phase plane for narrow band random excitation of a Dulling oscillator. J. Sound Vib. 104, 277 (1986).
508
A. H. NAYFEH and S. J. SERHAN APPENDIX
A
The matrix B and vector Q E
r
2c,
0
zc,
0
0
0
00
0
2c,
c3
c2
B=O
2c,
0
0
0
0
cz
0
0
0
0
E 2wo
c,+c,
0
” 00
0
B
c,
0 E
20,
2 00
0
 B
c,
c3
2
B
00
0
0
0
Cl 
00
0
0
0
c3
2
c2
B cq

2
(A3
.
APPENDIX
B
Method of equioalent linearization An equivalent
linear equation
that will furnish ii + o:u
an approximate
+ neti = e[hOcosRt
solution
 k,sinRt
The equivalent stiffness w.’ and damping pC can be determined difference between equations (1.1) and (Bl); that is o* = w2 + s(u& c 0
of equation
(1.1) can be written
+ c(t)].
as (Bl)
by minimizing
the meansquare
value of the
u)>  Fc. (B2) (u2>
(Iig(u,
Ii))  u$(uti)
P. = s
(B3)
where g(u, ti) =  ~5~ti + f~5~ti~ + au3. Next, we separate
u into a mean component
u. and a random
component
(B4) ut according
to
u=u,+u,
(BS)
where ic, + 0.’ u. + n&,
= s[hO cos Rt  k, sin fit]
ii, + cOful + /J&i, = Substituting
equation
(B4) into equations
Using equation
(B7)
(B2) and (B3) and using equation
~3.’ = 03; + [ui + (u:)]‘[&au: & = &cJ
+ [ti; + (#)I_’
(B6)
&S(t). + 6u*u~(u~)
(BS) yields + ra(ti:)]
[fac5,ti;: + 2&621i:(tii:)
+ fsS,(Ii;).
VW P9)
(B6), we obtain u; = +D’(h;
+ k;),
u; = $D4(h: + 2h;k;
ti; = R2u;
(BlO)
+ k;),
ti; = @u;
WI)
I”.
W2)
where D = [(co*c  S.l’)* + (p e fl)‘]
A shorttime averaging is implied in equations (BlO) and (Bl 1). For a Gaussian response (u:) = 3(u:>* and for a psuedo sinusoidal response (u:) = j(u:)‘. When < is a narrowband random excitation, the complex frequency response function H(o) can be obtained by combining equations (B7) and (3.27); that is H(w) =
EB”‘R
(0.’  co* + ip,o)(@
 0’ + i/?w)’
(B13)
Response statistics of nonlinear systems
509
The meansquare values of u, and J, can be written as
~2adi&.z + 87 + B(Q2 + r31
hi) =p.cof [(Q2  UJ:)* + 01. + 8) (J&n* + sol:,] (Ii?) =
~*~~&*(A+ B) r.C(Q2 4,’ + tP. + m.A~* + Bd)l~
0314)
PM
For the case of a white noise, we have (Uf)
=$, c
c
&*ns,
(Ii?) = .
A
(BlW