Journar
of Sound
THE
and Vibrarion
(1990) 139(l),
RESPONSE
FUNCTION
l-8
ENVELOPE
PROBABILITY
OF A DUFFING
RANDOM
DENSITY
OSCILLATOR
NARROW-BAND
WITH
EXCITATION
H. G. DAVIES ANI) Q. LIU Department
qf Mechanical
Engineering,
(Received
University of New Brunswick, Canada E3B 5A3
Fredericton,
New Brunswrck,
11 April 1988, and in revised .form 13 April 1989)
A Duffing oscillator excited by narrow-band random noise can exhibit multi-valued behaviour, with random jumps between two pseudo-stable states. In this paper an approximate form of the probability density function or the envelope of the response is obtained. State equations for the response envelope with the random excitation envelope as input are obtained by stochastic averaging. The joint probability density function of response envelope and input envelope can be found approximately for certain ranges of the system parameters. The response envelope probability density function shows two maxima in the cases where the response is multi-valued. Numerical simulation confirms the results.
1. INTRODUCTION
A Duffing oscillator between two distinct
excited by random narrow-band noise can exhibit random jumps slowly varying levels [l-4]. Equations that describe the time evolution have been obtained, and recently Rajan and Davies [4] have of the response envelope shown that these equations agree very well with numerical simulations, even in cases where multiple-valued response occurs. Although the envelope evolution equations are of interest in their own right, further information on the response process requires a knowledge of the statistics of the response, in particular the probability density function of the response. The probability density function for Duffing’s oscillator excited by white noise is well known, but no work appears to have been reported for narrow-band excitation of any type of non-linear oscillator, although Schmidt [5], for example, has presented approximate results for parametric excitation. An approximate probability density function is given here for the limiting case of very small excitation bandwidth. The approach is based on rewriting the envelope equations of Rajan and Davies [4], and solving, at least approximately, the associated Fokker-Planck equation.
2. EQUATIONS Consider
be modelled
a Duffing
oscillator
excited
OF MOTION
by a narrow-band
random
input.
The system can
by the equations j;+gIpj,+~~(y+FZ~y3)=&2f(t), f+E’r~+n;‘S=
The narrow-band excitation f is obtained f has centre frequency R, and bandwidth
E(rfl;)“W
(1) (2)
by filtering stationary Gaussian white noise W. T. The ordering of the parameters is such that 1
0022-460X/90/
100001
+
08 %03.00/O
@ 1990 Academic Press Limited
2
H. G.
DAVIES
AN13
Q.
LIU
E [.f’] = rr.9, and the spectrum off limits to rrS,,( 8(0 - 0, ) + a(0 -t f2, ))/2 as r tends to 0, so that the sinusoidal solution is obtained at this limit. E > 0 is a small gauge parameter that is required for the stochastic limit theorem to be applicable [6,7]. Equation (2) is just the linear form of the equation discussed by Roberts and Spanos [7]. The output of the filter f is used to drive the non-linear system. The system itself has centre frequency 0, and bandwidth /3, and a represents the strength of the non-linearity. Again the parameter E is introduced to order the terms. The case of most interest here is when 0, = 0,. The closeness of the excitation and system frequencies is emphasized by defining a detuning parameter: F26 = (0; - Ri)/(2L+).
(3)
As discussed in references [4] and [7], the responses f and y are expected to be approximately sinusoidal at frequency fif with random amplitudes and phases varying slowly with time scales much longer than the fundamental period 2a/R,. This feature is brought out by the transformations
The transformations
f=acos@,,
j‘=-&&sin@,,
@,=&r-+&J,,
(4)
y=bcos@,,
j = - 60, sin Q2,
@z=L?,-r+C$,.
(5)
(4) and (5) yield the following R,ci = - curare
exact equations:
sin’ @, - E( 70:)“2
sin @, W,
(6)
cos 0, W,
aR,&, = - e27-aR, sin @, cos @, - ~(70~)“~
(7)
0,b = - E2/3bflf sin’ CD,+ E2R&b3 cod Q2 sin Q2 - &*2bL?pScos Q2 sin a,-
&‘a cos @, sin Q2,
(8)
bR,&, = - &*pbL?, sin G2 cos Q2+ &*R&yb3 cos4 Q2 - E22bOf6 cos2 Gz - E’a cos @, cos Oz.
(9)
Equations (6)-(9) are in the standard form required by the Stratonovitch-Khasminskii limit theorem. Equations (6) and (7) are, of course, just the linear versions of equations given by Roberts and Spanos [7]. In the limit E + 0, (a, b, c$, , c$*) approaches a Markov process described by averaging equations (6)-(9) over a period 27r/L$ to yield the It8 equations, da=
E~{-(~(~/~)+(T~S~/~U)}~~-E(T~S~)”~~B,, db = &*{-(pb/2)+(a
de=e’{&+a
sin e/20,)}
dt,
(11)
cos 8/2R,b}dt-E{(T~S~)“~/.}~B,,
where f3= 4, - c$~, and B, and B2 are two independent The non-linearity is included in the term 0, = S -3f&b2/80,. For the as in the describing process b
(10)
unit-amplitude
(12) Wiener
process.
(13)
filter alone, the amplitude process c( is uncoupled from the remaining process, Roberts and Spanos case. The relative phase, 8, is clearly of importance in power input from the filter to the oscillator, so that the response amplitude cannot be similarly uncoupled.
RANDOM
RESPONSE
OF
A DUFFING
OSCILLATOR
3
3. COMPARISON WITH PREVIOUS WORK Equations similar to equations (lo)-(12) have been used before by Rajan and Davies [4]. With the definition M, = a2, mz= b’, P= a b sin 0, Q = ab cos 0, one obtains from equations (lo)-( 12) (noting that second derivative terms can be required in the transformation of It&type equations) dm, = E’(- Trn, +27r7So) dt - s(4m,nrSO)“’ dm, = E’(- Bmz+ P/0,)
dB,,
dt,
dP=e’(-(7+P)P/2+n,g+m,/(2nl))dt-&(~~SO/m,)”*(PdB,+QdB?), dQ=&Z(-(7+p)Q/2-~n,P)dt-&(7~S0/m,)”’(QdB,-PdB-).
( 14) (15) (16) (17)
In reference [4], the Wiener processes in equations (16) and (17) were not included. In the averaging process, Rajan and Davies [4] took only the mean and not the fluctuating component. They showed, however, that for small values of 7 the equations model the original equations (1) and (2) very well. Four equations were used even though, of course, the random variables are related by P2 + Q’ = m, m2. A “stationary” solution can also be obtained by taking expected values of mean square type terms in equations (14)-(17) and putting d/dt = 0. One then finds that 2B0f(0f(mz)+
?)m,=
?m,,
(18)
where ? = (~+/3)/2. This is a cubic equation for mean square response m, in terms of the mean square force m,. It gives the sinusoidal result of Nayfeh and Mook [8] in the limit T + 0 as well as the correct random linear result. It will be seen that equation ( 18) plays a role in the probability density function of b. Rajan and Davies showed that for very small values of T, the response follows the curve (18): i.e., an instantaneous value of m, can be calculated from the stationary equation, the response jumping if necessary from one stable branch of a multi-valued curve to another.
4. FOKKER-PLANCK EQUATION It has been found easier to obtain solutions of the approximate Fokker-Planck equation for the probability density function by first making yet another transformation of variables. Define x = - a cos 0, y = - a sin 0 and transform equations (lo)-( 12) from (b, a, @) to (b, x, y) variables. The new ItB equations are db = E’{-@b/2)
- (y/20,)}
dt,
(19)
dx=~2{-(~~/2)-~~y+(xy/2~nl-b)}dt-{~(~~S~)”’/a}(xdB,-ydB,),
(20)
dy=e’{-
(21)
(7y/2)+0;c-(x2/20,b)}dt-{&(T~S0)”2/a}(ydB,+xdB2).
Note that when the Fokker-Planck equation is formed, the result x’t y’ = a’ can be used to eliminate a, and gives an equation for p(b, x, y). In the linear case, that is with fl, a constant, the transformations r=y+flPR,b,
s=x-2pf&f&b/(T+p)
(22)
yield random variables (b, r, s) that are mutually statistically independent, so that the Fokker-Planck equation is easily solved. In the non-linear case, equation (22) describes a non-linear transformation. As might be expected, the variables (b, r, s) are no longer mutually independent. If an equivalent linearization of the non-linear term is made at this point, the resulting probability density function is Gaussian in r and s, and Rayleigh
4
H.
G.
IIAVIES
AND
Q.
LIU
in b. This form cannot model the bimodal character of the response when random jumps occur. An approximate solution is sought here by first deriving auxiliary equations from the Fokker-Planck equation, and postponing equivalent linearization until absolutely necessary. An approximate solution of the Fokker-Planck equation after just the linear transformation to r can be made by assuming that
p(b, r, xl =p,(b, x)pAr). The Fokker-Planck
equation
(23)
for p(b, r, x) is
-t($)+&[(+;K- ?r+R,x- x+~)p]-;(.,s,)($+$)=o. ixn,r+&+/3fl,*.b)p]
2Rfb
The substitution
(24)
(23) yields
pZ( r) = (const.) exp (- 7r2/ m&). One can also write p,(b, x) =p3(b)p4(x (b). The r” and r’ terms equated separately to zero to give
in equation
= - (2/ rrrSo)( TX - Pflri2eb)p4
dp,/dx
(24) are
(25)
and - ( P$f+P4$
+w-$(&-fk)P3P4}
-{
~ex-&+~}(~)p3p4=o.
(26) Equations (25) and (26) cannot be satisfied consistently in the non-linear case. The approximation is now made that the conditional probability satisfies, as it does in the linear case, the result ap,/ab This result
in equation
= (2pfl,-L?,/
TT&)(X - P12rR,b/ ?)p4.
(27)
(26) yields ap3/ab=(l/b)p~-(2P~~/~~So)(~nZ,+~’)bp3.
(28)
The assumptions are valid if equations (25), (27) and (28) can be satisfied. Clearly, equations (25) and (27) can only be satisfied in the linear case. Now use equivalent linearization to obtain an approximate solution by replacing bLL(b) in equations (25) and (27) by bfie, where 6, is a constant defined by ii, = E [b%,(b)]/E An approximate
solution
for p(b, r, x) can then be written
p(b, r, x) = (const.)p,( with the required written as
response
[b2].
b) exp {-(7/r&,)(x
envelope
p,(b)=(const.)bexp(
probability
(29) in the form
-/3L?,fii,b/f)* density
function
-~l(nt+i*)bdb).
- 7r2/mSo}, from equation
(30) (28) being
(31)
RANDOM
RESPONSE
OF
A DUFFING
OSCILLATOR
5
This form, of course, gives the correct Rayleigh value in the linear case (when approximations, inherent in the averaging approach, such as 0, = Q,, are made). Of more interest is the fact that this form can exhibit multiple maxima. The turning points of p,(b) are found by putting dp,/db = 0 from equation (31). They are at the roots of the equation 2p.n;b’(n;
+ ?) = +rSo.
(32)
This is essentially the same cubic equation for b2 as the “stationary” result given earlier. When the cubic has three real roots, corresponding to two stable and one unstable value of the purely sinusoidal response, the probability density function p,(b) exhibits two maxima and a minimum. If, for example, the bandwidth of the excitation is increased to a region where equation (32) has only one real root, then the probability density function p\(b) exhibits only one peak. Examples of these two extremes are shown in the numerical simulation results which follow.
5. NUMERICAL SIMULATION AND DISCUSSION Numerical techniques have been used to simulate the system described by equations one to approximate Gaussian process with zero mean and spectral density S(0) is a sum of cosine functions, (1) and (2). Among many available methods, a convenient
J(t)=(2a’/N)“J
; cos(R;r+B,),
(33)
i-l
5
0
5
b
Figure 1. (a) Numerical simulation of response envelope of the system given by equation (1) and (2) for e4S,, = 0.5, a = 0.5, p = 0.2, r = 0.02, 0, = 1.75, R, = 1.0 and N = 100. (b) Probability density function of the response envelope obtained by sampling the time history of Figure I(a). (c) Theoretical probability density function of the response envelope for the oscillator of Figure l(a).
6
H.
1.2.
,
,
,
(
,
,
G.
,
DAVIES
,
AND
Q.
LIU
, ib)
Figure 2. As Figure
-
1, but for ~=0.5
and N =400.
b Figure 3. Probability density function of response envelope C, T = 0.5. All other parameters as in Figure I(a).
for different
values of 7. A, 7 = 0.02;
B, 7 = 0.1;
where t$ are random angles with uniform distribution in (0,27r) and L$ are random frequencies distributed according to the density function g(0) = S(L?)/a’ with a2 = I”_: S(0) do. It has been shown [8] that W(t) tends to a Gaussian process when N tends to infinity. The procedure has been used a great deal and received detailed discussion in references [8,9]. For the case of concern here, expression (33) has been used to model the excitation f in equation (1). It is well known that the spectral density S(0) is described by s(n)=
E’~~~s,/{(n~-n*)*+(ET~)‘}.
(34)
RANDOM
RESPONSE
OF
A DLJFFING
OSCILLATOR
b Figure 4. Probability density function of response envelope C, S,, = 1.5. All other parameters as in Figure I(a).
for different values of S,,. A, S,, = 0.3; B, S,, = (1.5;
b Figure 5. Probability density function of response envelope C, a = 1.5. All other parameters as in Figure l(a).
for different
values of a. A, n = 0.1; B. a = 0.5;
Equation (1) with the excitation given by equations (33) and (34) can be integrated directly to give the time history y(r) and hence the response envelope b by squaring and using a low-pass filter. This direct simulation of the basic equations is shown in Figure 1 (a). The probability density function, p(b), is shown in Figures I(b) and I(c). Figure l(b) is obtained by sampling the time history shown in Figure l(a). It is understood that although a more accurate probability density function could be obtained by ensemble averaging, the cost of this was considered too high for use in this work. Finally, Figure l(c) is a plot of the analytic probability density function given in equation (31). The two maxima of the probability density function match the simulated results, although there are some differences in level. In Figures 2(a), (b) and (c) are shown results similar to those in Figures l(a), (b) and (c), respectively, but for a large value of 7, all other parameters being kept the same. It is clear that random jumps between levels do not occur, and the corresponding probability density function show only one maximum. Again agreement between simulation and analysis is reasonable.
8
H. G. DAVltS
ANI) 0. LIU
The effect on the probability density function of varying the parameters is shown in Figures 3-5. Figure 3 shows the effect of various bandwidths, Figure 4 of various excitation levels and Figure 5 of various non-linearities. In each case it is seen that either multiple maxima corresponding to multiple response levels and associated random jumps, or single-maximum behaviour can occur.
6. CONCLUSION
A Duffing oscillator excited by random narrow-band noise can exhibit random jumps between two pseudo-stable levels. A probability density function has been obtained for narrow bandwidth that shows two local maxima corresponding to the two pseudo-stable levels. The probability density function shows only a single local maximum for wider bandwidths.
ACKNOWLEDGMENT The work was supported of Canada.
by the National
Sciences
and Engineering
Research
Council
REFERENCES 1. T. KAPITANIAK 1985 Journal of Sound and Vibration 102, 440-441. Stochastic response with bifurcations to nonlinear Duffing’s oscillator. 2. H. G. DAVIES and D. NANDLALL 1986 Journal ofSound and Vibration 104, 277-283. Phase plane for narrow band random excitation of a Duffing oscillator. 3. T. FANG and E. H. DOWELL 1987 International Journal of Nonlinear Mechanics 22, 267-274. Numerical simulations of jump phenomena in stable Duffing systems. 4. S. RAJAN and H. G. DAVIES 1988 Journal ofSound and Vibration 123,497-506. Multiple time scaling of the response of a Duffing oscillator to narrow band random excitation. 5. G. SCHMIDT 1976 Proceedings of the VII International Conference on Nonlinear Oscillations, Berlin, 341-359. Nonlinear systems under random and periodic parametric excitation. 6. R. L. STRATONOVICH 1963 Topics in the Theory of Random Noise. New York: Gordon and Breach. 7. J. B. ROBERTS and P. D. SPANOS 1986 International Journal of Nonlinear Mechanics 21, 111-134. Stochastic averaging: an approximate method of solving random vibration problems. 8. M. SHINOZUKA 1971 Journal of the Acoustical Society of America 49, 357-367. Simulation of multivariate and multidimensional random processes. 9. M. SHINOZUKA 1972 Computers and Structures 2, 855-874. Monte Carlo solution of structural dynamics. 10. L. ARNOLD 1974 Stochastic Di$eren?ial Equations. New York: John Wiley.