- Email: [email protected]
Transient response of non-linear systems to random excitation
Journal of Sound and Vibration (1981) 74(l), 1l-29
TRANSIENT
RESPONSE
OF NON-LINEAR
TO RANDOM
SYSTEMS
EXCITATION
J. B. ROBERTS School of Engineering
and Applied Sciences, University of Sussex, Falmer, Brighton BN19QT, England (Received 6 June 1980)
An approximate method for computing the transient response of non-linear oscillators to suddenly applied, wide band random excitation is developed. The technique is derived initially for systems with non-linear damping, by using the method of stochastic averaging, and then generalized to deal with oscillators with arbitrary non-linear restoring forces. Predictions from the theory are compared with corresponding digital simulation estimates, in a number of specific cases.
1. INTRODUCTION
In many random vibration problems it is reasonable to model the excitation and response processes as stationary random processes. However, it is important to realize that such modelling is always an approximation, since all real excitation processes are nonstationary to some extent. In some applications there is a very severe degree of nonstationarity in the statistical character of the input which may be best represented as an abrupt change. Examples of such non-stationary excitations are wind gusts, earthquakes, ignition transients in missiles and spacecraft, and aircraft landing loads. In all such cases one is concerned with studying the transient response of structural systems to stochastic excitation. The simplest transient problem, and that which has been most extensively studied, is to compute the statistics of the response of a structural system to suddenly applied, stationary excitation, assuming that the system is initially at rest. For linear systems this problem can be solved fairly easily by a variety of methods, most of which are based on the convolution and Fourier transform relationships (e.g., see references [l-13]). A number of timedependent statistics, such as the mean square response us. time, have been computed by these methods for single degree of freedom systems and for some multi-degree of freedom systems. All real systems are to some extent non-linear, however, and it is thus a matter of practical concern to develop methods of computing the transient response of non-linear systems to stochastic excitation. Unfortunately, it is difficult to extend the standard linear methods to the non-linear case. For the computation of the mean square response Sakata and Kimura have shown that it is possible to extend their method [ 131, for non-stationary excitations, by using the method of equivalent linearization (e.g., see reference [14]), although there are very considerable complications in this approach [ 151. In principle the perturbation method (e.g., see reference [16]) can also be extended to the non-stationary case but again the details are complicated and such an approach has not yet been attempted. Faced with these difficulties, several workers have resorted to the direct digital 11 0022-460X/81/010011
+ 19 $02.00/O
@ 1981 Academic Press Inc. (London) Limited
12
J. B. ROBERTS
simulation approach to obtain estimates of non-linear system response (e.g., see references [17] and [18]). However, when the stochastic excitation is wide band in character an alternative theoretical approach is possible, which is not derived from linear methods. This involves setting up an appropriate Fokker-Planck Kolmogorov (FPK) diffusion equation for the transition density function, and solving this to find the transient response. Unfortunately, even for single degree of freedom systems, the transient response cannot be found analytically and a numerical solution procedure must be used. This approach has been used by Wen [19] to obtain various response statistics for the Duffing oscillator. An alternative numerical procedure has been given by Toland and Yang [20] and Toland et al. [21]. In the FPK approach it is possible to simplify the analysis considerably if the system is a lightly damped oscillator with a small degree of non-linearity. One can then exploit the fact that the amplitude envelope process is approximately a Markov process, governed by a one-dimensional FPK equation (e.g., see reference [22]). This one-dimensional equation is much easier to solve than the exact, two-dimensional FPK equation [19-211. In particular, one can readily find the time-dependent transient response of the amplitude envelope, to a non-stationary wide band input, from the FPK equation. Such results have been obtained for a variety of non-linear oscillators by Iwan and Spanos [23, 241 and Spanos [25-291, using a perturbation technique to find an approximate series solution in terms of eigenvalues and eigenfunctions. It is the purpose of this paper to show that the FPK method, as applied to the amplitude process of a non-linear oscillator, yields considerably more information than just the statistics of the amplitude process itself. Firstly it is demonstrated that it also yields an approximation to the time-dependent joint distribution, p(x, i ; t), of the displacement and velocity processes: from this distribution a wide variety of statistics can be calculated, such as level crossing probabilities. This part of the paper is essentially a generalization of the analysis given in reference [22]. Secondly, it is shown that this approach can be further generalized to deal with oscillators with arbitrary non-linear restoring forces, by using the fact that the energy envelope is Markov [30-331. By solving the FPK equation for the energy envelope it is shown that one can obtain not only the statistics of the energy envelope but also an approximation to the joint distribution p(x, f ; t). For the special case of oscillators with restoring forces governed by a power-law relationship, a simple analytical solution is obtained, which generalizes the previously known result for a linear oscillator. In other cases the appropriate FPK equations are solved by using a random walk analogue [34,32] to form a simple numerical scheme. To assess the accuracy of the approximate theory which is developed here, results obtained from it are compared with corresponding digital simulation estimates, for a variety of cases.
2. THE LINEAR
OSCILLATOR
To introduce the proposed methods it is convenient to begin with the case of a linear oscillator, since here exact solutions are obtainable by standard procedures. The appropriate equation of motion is i+225ooi
+&x
= y(i),
(1)
where n(t) is the displacement response, l is the non-dimensional damping factor and o. is the undamped natural frequency. y(t) will be assumed to be a stationary random process, with zero mean and of wide-band character, which is suddenly applied (i.e., “switched on”) at time t = 0.
TRANSIENT
NON-LINEAR
SYSTEM
13
RESPONSE
The spectrum of y(t) may be defined as
S(o)=& 1-1 W(T)ei‘“’d7,
(2)
where w(7) is the covariance function for y(t): i.e., (3)
w(~)=E{YWY(~+~)I.
It will be assumed that S(o) is smooth in the neighbourhood of w = 00 and that the damping is light. In these circumstances y(t) can be well approximated as a white noise process, n(t), with spectrum S,(w) = S(wo) a constant [l]. The covariance function for n(t) is then W”(7) = 2&(oo)S(7). (4) With y(t) modelled as a white noise it is fairly easy to compute the mean square response of the linear oscillator, CT:,by using the convolution integral and impulse response function (e.g., see reference [l]). The result is
&> -=
1-e-2w
cl2
l+_ [
x2
l-l2
.
2
sm
l@O
_
u’+dl-gz
-sin2fGt
1
(5)
, t>O,
where 6 = 6JoJl
(6)
is the damped natural frequency and CT’is the stationary mean square response, given by Cr2= 7rs(wo)/256J:.
(7)
As 5 -* 0, equation (5) approaches the asymptotic result &(t)/a’
= 1 - e-2yoo’.
(8) In Figures l(a) and (b) r:(t)/ a2 is plotted against the non-dimensional time CN, where N = oot/2r (the number of “cycles”), for [ = 0.05 and O-20, respectively, according to I
Statlanary
(b)
response
Statlonory
I
I
response
0.7 -Jb
0.6
N‘ b’
05 0.4 0.3 0.2 0.1
0
0.1
02
0.3
0
0.1
0.2
O-3
CN Figure 1. Variation of the mean square displacement with time for the linear oscillator. Comparison between the exact solution and the asymptotic solution. (a) [ = O-05; (b) 5 = 0.20.
14
J. B. ROBERTS
equation (5). In both cases a comparison is made with the asymptotic result of equation (8). It is seen that the effect of the sine terms in equation (5) is to produce an oscillation, of frequency 26, superimposed on the smoother variation given by the asymptotic result. The amplitude of this oscillation is of order 5 and the oscillation is seen to be fairly insignificant in the case [ = 0.05. It can be concluded that, for light damping (5 < 0.1, say) of the kind usually encountered in engineering structures, the asymptotic result is a good approximation. This suggests that, in developing a suitable approximate method for lightly damped oscillators, which may not necessarily be linear, one should apply some kind of averaging, or smoothing, procedure, which will eliminate relatively small oscillatory terms such as those appearing in equation (5). Such a “stochastic averaging” procedure has been developed by Stratonovitch [35], and is a generalization of the well-known deterministic averaging method of Bogoliubov and Mitropolski [36]. The averaging method is based on the amplitude envelope process, A(r), where A*(t) = x*(t)+i*(t)/&
(9)
which has a simpler statistical structure than either x(t), or i(t). This can be demonstrated quite easily by considering the mean square amplitude E{A*(t)}. Using the impulse response function again one finds that the mean square velocity, a,*, is given by
af0)
a20Zo=
l-e-2co~r
[
x2
l+-
.
l-l2
2
s1n
iSJO
_
-sin2Gt
“Q-q2
.
I
t>O.
(10)
Hence, from equation (9), w*w1=
I_
1 +
e-2&o,r
2a2
[
x2
sinZ
1_52
&
I *
(11)
The oscillatory term here is of order t*, which will thus be extremely small when the damping is light (cf. equations (5) and (lo), where the dominant oscillatory term is of order 4’). Thus a very good approximation to equation (11) is E{A*(t)} = 2g*(l -e-2rWot)
(12)
and one thus has the asymptotic result (see also reference [25]) a: = &&
= i~{A*(t)}.
(13)
This shows, in a special case, how the smoothed (or averaged) behaviour of the statistics of x(t) and i(t) are related to the statistics of A(t). In the next section a more general relationship between x(t), i(t) and A(t) will be derived, for a class of non-linear oscillators which includes the linear case. 3. THE STOCHASTIC
AVERAGING
METHOD
Consider now a non-linear oscillator, with an equation of motion of the form i+@F(x,k)+&x=y(t),
(14)
where /3 is a coefficient which is small in magnitude. Again y(t) is assumed to be a wide band process, with zero mean. Since the energy dissipation will be small (@small) x(t) and i(t) will have an oscillatory, or narrow band, character and it is meaningful to define an amplitude process, A(t), and a
TRANSIENT
NON-LINEAR
SYSTEM
RESPONSE
15
phase process, 4(t), as follows: x=Acos@,
li=--woAsin@,
(15)
where @=oot+q!L
(16)
A(t), x(t) and i(t)
are then related according to equation (9). The equation of motion can be recast as two simultaneous equations in A(t) and 4(t). On applying the averaging procedure due to Stratonovitch [35] to these equations (see the Appendix in reference [22]) one eventually obtains the approximate equations A = -(p/wo)G(A)
+{~S(wo)/2Ao~}-[nS(wo)]“‘n(t)/wo,
4 = -(PIAwo)H(A)
(17)
- W(oo)11’2N)lAoo,
(18)
F(A cos @, -w. A sin 0) sin ad@,
(19)
F(A cos @, -ooA
(20)
where 2.77
G(A)=-l
2lr I 0
and n(t) and f(t) are mutually independent, such that E{n (t)n 0 + 7)) = WW(t
sin @) cos @ d@,
stationary Gaussian white noise processes,
+ r)] = S(T),
E{n(t)&+T)}=O.
(21722)
It follows from equations (17) and (18) that the joint process A(t), 4(t) is a twodimensional Markov process with a transition density function p(a, c$]u~,c$~;t) governed by the FPK equation
aP a at au
-&-
PWU) ap+dbo) Mwo) %4-y--~P +II au0 a4 -2+7y. Uw. 20~ 0
a2p 1 a2p [ aa
a
a4
3
(23)
(For the special case H(A) = 0, this agrees with the result given in reference [22].) Furthermore, it follows from equation (17) that A(t) is a one-dimensional Markov process, with a transition density function, p(u]uo; t) governed by the FPK equation
ap a -=at
aa
&?(u)-$$$}p]+%$ o. 0
(24)
Here the transient response is sought, with the oscillator assumed to be initially at rest, at time t = 0. The time dependent joint distribution of A(t), 4(t) is then given by
0, 4; t) =p(a, 41% 40; 0,
(25)
where p must satisfy equation (23). On comparing equations (23) and (24) it is found that the axisymmetric form Pb,
4; t)=pb;
tmr,
(26)
where p(u; t) = p(u (0; t) satisfies equation (23). Thus one has the remarkable result that the time dependent solution of equation (24), p(uIO;t), also gives the solution to equation (23). This only applies, of course, with the initial condition A0 = 0.
16
J. B. ROBERTS
From equation (26), by transformation of variables from a, 4 to x, i’, one can express the joint distribution p(x, i; t) of x(t) and i(t) in terms of ~(a; t). The result is P(& i; t) = (claMa;
t),
(27)
where a = (X2+.P/o~)1’2
(28)
and c is a normalization constant. Equation (27) can be used to compute various statistics of x(t) and i(t). For example, the mean square of r(t) is given by -x
dxdi=c
u*p(u ; t) cos* @ dud@ = E{A*(t)}/2,
(2% which is in agreement with equation (13). 3.1. SOLUTION OF THE FPK EQUATION To find ~(a; t), and hence p(x, i; t), the time-dependent required, subject to the initial condition
solution of equation (24) is
p(u; 0) =6(u).
(30)
Unfortunately, a closed form analytical solution seems possible only in the linear case. This well known result is as follows (e.g., see reference [35]): a pb;d=
2
wo(1-4)
exp
(
-
a
2
2&(1-4)
(31)
I
where 4 = exp
cu=a/a
(-25w0tL
W33)
and CTis given by equation (7). Hence, from equation (29), it is readily found that af(t)/a*
= 1 - e-2Loof,
which, as previously shown, is asymptotically correct as [+ 0 (equation equation (31) represents simply a time-dependent Rayleigh distribution, ~(a’; t) = a’ exp (-a’*/2)
(34) (8)). Hence, (35)
where (Y’= a/&&),
(36)
as Spanos has noted [26]. For the non-linear case a numerical solution of equation (24) is necessary. For this purpose the random walk analogue discussed in reference [34] is useful. This involves replacing A(t) by a discrete random walk process, R(ti), which can only assume the amplitudes & = k Su (k = 0, 1, . . .) at time tj = j 8tG = 0, 1, . . .). This process is such that, from state (&-, it can only move to either state (&+I, ti+l) with probability pk, or to state (a&l, ti+l) with probability qk = 1 -pk. With an appropriate choice of pk it can be shown that R(tj) is governed by the same diffusion equation as A(t) (equation (24)), in the limit as St + 0, provided that tj)
@a)* = {7rS(oo)/o~} (at).
(37)
TRANSIENT
NON-LINEAR
17
SYSTEM RESPONSE
In practice the theoretically infinite range of ak must be truncated by choosing some upper limit to k (M say); it is then convenient to set an absorbing barrier at the level UM. The computed distribution, p(& ; ti) will be virtually independent of the value of M, if this is sufficiently large. At ak = 0 the random walk process becomes singular, but this difficulty can be overcome by placing a reflecting barrier at this level. Results obtained from this numerical method will be compared later with corresponding digital simulation results.
4. OSCILLATORS
WITH
ARBITRARY
NON-LINEAR
RESTORING
FORCES
The stochastic averaging method, as described above, is not a useful technique for examining the effect of non-linear restoring forces on the transient response. For example, consider the case of a Duffing oscillator with linear damping; the F(x, i) function in equation (14) is then F(x, i) = i + E&X3.
(38)
On evaluating the G(A) function (see equation (19)) it is found that the cubic term does not contribute, i.e., (39)
G(A) = ooN2,
as for the linear oscillator. Iwan and Spanos have attempted to overcome this difficulty by combining an equivalent linearization method with stochastic averaging, and treating the equivalent frequency as amplitude dependent [23,25,26]. However, they include terms which are not consistent with the order of approximation inherent in the stochastic averaging method and the validity of their approach is doubtful [37]. Here an alternative approach is followed which is based on a consideration of the energy envelope. Consider an oscillator governed by the following equation of motion: f+/?F(x,1)+G(x)=y(t), where G(x) represents an arbitrary process, V(t), is defined as
non-linear
V(r) = P/2 + U(x),
(40)
restoring force. The energy envelope
U(x) = ix G(t) dS.
(4~42)
Jo
Here i2/2 is the kinetic energy of the oscillator and U(x) is its potential energy. In the case where G(x) =&x, previously considered, U(x) = o&*/2 and then V(t) = &A*(t)/2, where A(t) is defined by equation (9). It has been shown previously by the author that V(t) is approximately a one-dimensional Markov process when p is small [30-331. The appropriate FRK equation for the transition density function p(u]uo; t) of V(t) is
ap a -$= ~[BBbd-
d%JoxlP~+
Moo)
gz
(43)
[C(U)Pl,
where 1 B(u) = &A(u)
F(x, J2(v - U)) dx,
Jv-vdx,
(44)
18
J. B. ROBERTS
(45) The integration range, R, in these integrals is such that V(x) < ~1. By solving equation (43), with the initial condition o. = 0, one obtains the timedependent distribution of V(t), p(v; t) =p(vlO; t),
(46)
with the initial form p(u; 0) = S(v). This gives useful information on the rate at which stationarity is achieved. However, one can go considerably further than this by using the fact that, for small p, the joint density of x(t) and V(t), p(x, u; t), is related to p(u; t) according to [30] p(x, u; t) =p(v; t)/2A(o)&
- U.
Upon reverting to x, f variables, equation (47) can be readily transformed expression for the joint density of x(l) and i(t), i.e., p(x, f; t) =p(v;
t)/hA(u).
(47) into an (48)
In the case of linear damping equation (48) agrees with the known exact result for the stationary distribution. Furthermore, in the case where G(x) = wfx then p(u; r) =P(a;
r)l&,
(49)
where tr = w&‘/2, and equation (48) reduces to equation (27). Although equation (27) was derived by using the stochastic averaging procedure, which has a rigorous foundation [38], a correspondingly rigorous derivation of the generalization given by equation (48) seems difficult to achieve. The derivation given in reference [30], whilst physically plausible, is admittedly rather heuristic in nature. A rather obscure argument in support of equation (48) has also been given by Stratonovitch [35]. Some physical insight into the validity of equation (48) can be obtained by considering the phase plane: x, f. In the absence of excitation, and with p = 0 in equation (40), trajectories in the phase plane will be simply orbits corresponding to a fixed value of u =1*/2+ V(x). When the energy dissipation per cycle is small (p small) and the excitation is correspondingly weak, each “cycle” in the phase plane will approximate (in a macroscopic sense) to a constant ~1orbit. For a suddenly applied input the orbit level will gradually increase from zero. In these circumstances it is reasonable to expect that a constant D contour in the phase plane is also a contour of constant probability density, i.e., P(X, i; r) =f(u; 0,
(50)
wheref(v; t) is some function of 0 and t. It follows from equation (50) that the joint density of V(t) and x(t) is given by p(u, x; r) =f(v; t)/_.
(51)
On integrating with respect to x one has p(u;t)=
OJp(v, x; t) dx =f(o; t)&A(v); I -al
(52)
hence (53) Combining equations (50) and (53) one recovers equation (48).
TRANSIENT
NON-LINEAR
SYSTEM
19
RESPONSE
By integration of equation (48) one can find various statistics associated with x(t) and x(t). For example, consider the evaluation of the mth moment of x(t), E{x”(t)}, E{x”(t)}=
Irn Irn x”p(x, x; t) dx dl --g, -ca
x”p(x, u; t) dx dv
=
where
4.1.
ANALYTICAL
SOLUTION
The FF’K equation for V(t) (equation (43)) can be solved analytically in the case where the damping is linear (i.e., F(x, II) = x) and the restoring force has the power-law form
G(x) = k]x]” sgn (x).
(56)
For this purpose it is convenient to define the non-dimensional
energy variable
Z(t) = v(t)/yka”+’
(57)
where (58) and (T is the stationary standard deviation of x(t). Then, as shown in reference [31], the solution to equation (43), when expressed in terms of the transition density p(z]zo; t) of Z(t), is as follows:
l p(zlzo;t>= (l_q)($Jp'(*)exP
{-‘;;_yp’ 51, )Ip(
(59)
where P = (1/a) - 1,
q ; e-olpr
(Y=2(v+l)/(v+3),
(60)
and I,( ) is the modified Bessel function. Here p(z; t) only is required; hence, setting z. in equation (59) to zero and using the fact that I,(x)~{l/~(l+p)}(x/2)p
asx+O
(61)
it is found that p(z; t> =
1 r(1 +p)(l -q)‘+‘zp
-_ f
exp (
l-q > -
(62)
In the linear case p = 0, (Y= y = 1, /3 = 24’~~ and equation (62) is equivalent to equation (31). It is interesting to note that the shape of the distribution of Z(t), as given by equation (62), is invariant with respect to time. Thus, if z’ = z/(1 -4)
(4
then equation (62) can be written as p(z'; t)= {i/r(i+p)}(~')p
exp
(-z').
(64)
20
J.B. ROBERTS
In the special case where t + co, then q + 0, z'= z,and equation (64) gives the stationary density function. This conclusion is a generalization of that given earlier for the particular case of a linear oscillator (see equations (31), (35) and (36)). The invariance with respect to time of the form of the distribution of energy, found here for the case of the power law oscillator, is due to the fact that B(n) and C(U), in equation (43), are identical and proportional to U. It is shown in the Appendix that time-invariance of the shape of p(v; t) occurs only when B(u) and C(v) are both proportional to ZI. The moments of Z(t) can be easily computed from the above results by suitable integration. It is found that E{Z”(t)}=(l-q)(p+m)(p+m-1).
. . (p+l),
(65)
and in the special case where m = 1 this reduces simply to E{Z(t)) = (1 -q)(p + 1) =
(1--4)/a.
(66)
In Figure 2(a) the non-dimensional mean energy E{V(t)}/k~‘+” = yE{Z(t)} is plotted against the non-dimensional time, pt, for Y = 0, 1, 2 and 3, according to equation (66).
Figure 2. (a) Variation of the mean energy with time for the power-law oscillator; effect of power index; (b) variation of the mean square displacement with time for the power-law oscillator; effect of power index.
Similarly the moments of x(t) can be found by combining equations (54) and (62), and integrating. For example, the mean square response of x(t) is found to be given by the following remarkably simple formula: &(t) = a*(1 - q)2'(1+v).
(67)
For the linear case (Y = 1) this result reduces to equation (34). Figure 2(b) shows the variation of &(t)/a* with fit for y = 0, 1, 2 and 3, according to equation (67). 4.2. NUMERICAL SOLUTION OF THE FPK EQUATION In general equation (43) cannot be solved analytically and a numerical solution procedure is required. For the results to be presented in the next section an extension of
TRANSIENT
NON-LINEAR
SYSTEM
RESPONSE
21
the random walk analogue has been used, as described in reference [32]. Briefly, the method is based on a non-linear co-ordinate transformation, which transforms equation (43) into a diffusion equation of the same form as a limiting random walk (at + 0). The co-ordinate transformation must be performed numerically but need be carried out only once, prior to the numerical diffusion process. The time step and amplitude increment must be related, in a similar way to that indicated in equation (37).
5. COMPARJSON
WITH
SIMULATION
RESULTS
To test the approximate theory developed a number of comparisons have been made between theoretical predictions and simulation estimates, as follows. 5.1,
SIMULATION
METHOD
The simulation method consists of generating a number of realizations of the transient response, from t = 0 to t = T, by digital means. For this purpose the input process y(t) is approximated as a step function which jumps from one level to another at equi-spaced time instants, at spacing dt. The levels are proportional to a series of independent, Gaussian random numbers. Such a step process is a very close approximation to a white noise process when At is small. The equation of motion is integrated by using a standard fourth order Runge-Kutta algorithm, with an integration step length of At. To obtain reasonably stable estimates of parameters such as the mean energy and mean square displacement, 4000 realizations were generated, for every case. These realizations were ensemble averaged to provide the required statistical estimates. 5.2.
OSCILLATORS
WITH
NON-LINEAR
DAMPING
It was shown earlier that, for the case of linear damping, the stochastic averaging theory gives a result for the mean energy envelope which has an error of order [*, whereas for the mean square displacement the error is of order f: One can thus anticipate that, in the more general, non-linear case the approximate theory is likely to give more accurate estimates of total energy statistics than of displacement (or velocity) statistics. To test the theory for non-linear damping an oscillator with a cubic damping is considered-the equation of motion is as follows: ~+2~~o(~+E1~3)+W~X=y(t).
(68)
If u. is the standard deviation of the stationary response in the linear case (&r = 0) then equation (68) can be conveniently non-dimensionalized in terms of X=x/u,
and
r =oot.
(69)
It then becomes X+25(X+eEX3)+X=
Y(r),
(70)
Y(r) = Ywbo&.
(71)
where &T= “;I_&, The corresponding
non-dimensional
energy is then v* = x2/2 +X2/2.
(72)
In Figure 3(a) the theoretical variation of E{ V*} with non-dimensional time cN(N = oot/27r) is shown, for ET =0, O-2 and 1.0. These result? were obtained by using the random walk method to solve equation (24) numerically. Also shown in this figure are
22
J. B. ROBERTS
0
1
I
I
/
O-04
0.06
042
M6
I
SN
Stationary ---------------
response
(
0.6 05 0.4 0.3 0.2 0.1
0
I
-+-
I
4
N Figure 3. (a) Variation of the mean energy with time for an oscillator with a cubic damping term, comparison of theory with simulation results; simulation estimates: [ = 0.04: v, ET = 0; X, ET = 0.2; +, ET = 1.0; (= 0.10: 0, ET = 0; 0, B: = 0.2; 0, ET = 1.0; (b) variation of the mean square displacement with time for an oscillator with a cubic damping term, C = 0.04, comparison of theory with simulation results; theory, -, simulation; (c) variation of the mean square displacement with time for an oscillator with a cubic damping term, 5 = 0.10, comparison of theory with simulation results; -, theory, m, simulation.
corresponding simulation estimates of E{ V*}, for 4’= 0.04 and 4’= O-10. The time step, AN = woAt/27r, in the simulation was 0.02. There is good agreement between the theoretical curves and the simulation estimates at all three values of ET. Figure 3(b) shows the computed variation of the mean square of X(r), (T;, with N for 4’= 0.04. It is noted that these theoretical curves can be obtained directly from Figure 3(a) since, according to the stochastic averaging theory, r&(f) = E{ V*} (see equation (29)). The corresponding simulation results vary in an oscillatory way with time, as one might expect in view of the exact linear result (see Figures l(a) and (b)): to bring out this characteristic the
TRANSIENT
NON-LINEAR
SYSTEM
RESPONSE
23
simulation results are here plotted as a continuous line, rather than as a series of points, as in Figure 3(a). Again it is noted that good agreement is obtained between the two sets of results. One would expect the exact solution for ET > 0 to be oscillatory, as in the linear case, and this is indicated by the simulation results, although this behaviour is to some extent obscured by the inherent statistical uncertainty. A similar comparison is made in Figure 3(c) for the case 5 = 0.10. There is a tendency for the amplitude of the oscillations in the simulation results to be greater at the higher damping factor: again one can anticipate this by reference to the linear case and also from the nature of the stochastic averaging theory. 5.3.
OSCILLATORS
WITH
NON-LINEAR
STIFFNESS
5.3.1. Power-law oscillators with linear damping Again it is convenient at the outset to non-dimensionalize is
the equation of motion, which
i+pi+k(xj”sgn(x)=y(t).
(73)
If X=x/u
and
(74)
r’=Pt
then equation (73) can be written as X+X If a non-dimensional
+ (k~“-‘/~2)~X~” sgn (X) = y(t)/&.
(75)
damping coefficient is defined as x = /3/2[ka”-‘1”’
(76)
(in the linear case x = l), then it is evident that the non-dimensional response depends only on Y and x. Figure 4(a) shows the variation of the mean non-dimensional energy yE{Z(t)) with time for Y = 0 and v = 1 (see also Figure 2(a)). These were computed from equation (66). Also shown are simulation results for ,y = 0*04,0*10 and O-50. At the highest of these damping values there is some discrepancy between the simulation results and the approximate theoretical results, but the agreement is much better at the lower values of damping. Figures 4(b) and (c) show comparisons between theory and simulation estimates for the mean square displacement, a& us. time, for x = 0.04 and x = O-10, respectively. It is clear that the agreement is best at the lower value of damping as one would expect. The oscillations in the simulation results are only partly due to statistical variability: there is little doubt that the exact solution would be oscillatory, the ripple having a frequency and amplitude similar to that exhibited by the simulation results. 5.3.2. Oscillators of the Dufing type Finally an oscillator of the Duffing type has been considered; i.e., one with an equation of motion as follows: i +2&k Again, non-dimensionalization one obtains
+0&L +&2x3) = y(t).
in terms of X and T, as given by equations (69) is useful: X+21X+X+&X3=
where Y(7) = y(t)/ao&
(77)
ET = &$EZand differentiation
Y(t),
is now with respect to T.
(78)
24
J. B. ROBERTS
I.0
0.9
0.8 0.7 0.6
N
2
0.5 0.4 0.3 0.2 01
0
0.1
0.2
03
0
0.2
03
04
P, z
Figure 4. (a) Variation of the mean energy with time for the power-law oscillator, comparison of theory with simulation results; simulation estimates: Y = 0: ?? , x = 0.04; +, x = 0.10; A, * = 0.50; Y = 3:0, x = 0.04; X, x = O-10; V, x = 0.50; - - -, stationary response; (b) variation of the mean square displacement with time for the power-law oscillator ,y = 0.04, comparison of theory with simulation results; (c) variation of the mean square displacement with time for the power-law oscillator, x = 0.10, comparison of theory with simulation results.
In Figure 5(a) the variation of the mean of the total non-dimensional v* = R2/2 +x2/2
+ E:x4/4,
energy, (79)
with time is plotted for E; = 0,0*2 and 1.0. These were computed from the FPK equation for V(t) (equation (43)), by using the extended random walk analogue method. Simulation estimates for E{ V*}, for 5 = 0.04 and O-10, are seen to be in good agreement with the theoretical curves. It is interesting to note that the total energy is only weakly dependent on the extent of the non-linearity. Figures 5(b) and (c) show comparisons between theory and simulation estimates for the mean square displacement, aat), for [= 0.04 and [ = 0.10: respectively. Again the
TRANSIENT
(b) I.0
NON-LINEAR
I
I
I
I
I
2
3
4
25
SYSTEM RESPONSE
_-___-__---___--__----------
0.9
0.8 D7 0.6 “2
0.5 0.4 0.3 0.2 @I
, 0
I I
I
I
I
2
3
4
N
Figure 5. (a) Variation of the mean energy with time for the Duffing oscillator, comparison of theory with simulation results; simulation estimates: 5 = 0.04: 0, ef = 0.2; 0, ei: = 1.0; 5 = O*lO: X, e: = 0.2; +, E; = 1.0; (b) variation of the mean square displacement with time for the Dulling oscillator, C = 0.04, comparison of theory with simulation results; - - -, stationary response; -, theory; MA,simulation; (c) variation of the mean square displacement with time for the Duffing oscillator, c = 0.10, comparison of theory with simulation results. - - -, stationary response; -, theory; m, simulation.
simulation
results
indicate
an oscillatory
solution
which
is close to the present
theory
when
the damping is light. 5.4. ACCURACY OF THE COMPUTED RESULTS It is interesting to enquire briefly into the accuracy of the computed theoretical results, obtained by using the random walk analogue. For this purpose it is convenient to examine some results obtained for the Duffing oscillator; for this case the stationary mean square
26
J. B. ROBERTS
AX Figure 6. Variation of the percentage error in the mean square stationary displacement, as computed by the random walk analogue, with the amplitude increment.
displacement can be computed exactly and compared with the asymptotic result obtained from the numerical diffusion process, as the elapsed time becomes large. In the extended random walk method the range for X(t) is,.initially, divided up into a number of equi-spaced amplitude increments, dx say. Thus dx is a measure of the coarseness of the discretization in the random walk method. Figure 6 shows the variation of the percentage error in the stationary solution, with dx, for Ez = 0.20 and E: = 1.0. This shows that, in both cases, there is an almost linear variation of error with Ax. In the results given earlier a value of Ax = O-01 was chosen; it is clear that the error is then only a few percent, at most. As a final point, it is noted that one can use the known, exact result for the stationary response to correct the results obtained by the numerical diffusion procedure. Of course, I
06
I
I
I
I
-
_.-.------. 0.5 -
04N
c
0
02
c-4
06
k0
Figure 7. Variation of the mean square displacement with time for the Dutling oscillator. ET = 1.0. Effect of correcting the results from the random walk method by using the known exact result for the stationary case. corrected estimate, Ax = 0.01. -* -, Original estimate, Ax = O-07; - - -, corrected estimate Ax = 0.07; -,
TRANSIENT NON-LINEAR SYSTEM RESPONSE
27
for the transient response such a correction will be an approximation. Figure 7 shows the computed variation of crgwith [N for ~5 = 1.0 and an amplitude increment Ax = O-07. For the stationary response the error is 19.6% (see Figure 6). If the ordinates are reduced everywhere by 19.6% one obtains a curve which is in very close agreement with the corrected result for Ax = O-01 (error 2.8%).
REFERENCES 1. T. K. CAUGHEY and H. J. STUMPF 1961 Journal ofApplied Mechanics, American Society of Mechanical Engineers 28, 563-566. Transient response of a dynamic system under random excitation. 2. R. L. BARNOSKI and J. R. MAURER 1969 Journal of Applied Mechanics, American Society of Mechanical Engineers 36, 221-227. Mean-square response of simple mechanical systems to nonstationary random excitation. 3. R. L. BARNOSKI and J. R. MAURER 1973 Journal of Applied Mechanics, American Society of Mechanical Engineers 40, 73-77. Transient characteristics of simple systems to modulated random noise. 4. J. B. ROBERTS 1971 Journal of Sound and Vibration 14,385-400. The covariance response of linear systems to non-stationary random excitation. The covariance response of 5. J. B. ROBERTS 1971 Journal of Sound and Vibration 17,299-307. linear systems to locally stationary random excitation. p. T. K. HASSELMAN 1970 Journal of Applied Mechanics, American Society of Mechanical 189. A comparison of solutions for the mean square response of a simple Engineers 37,1187-l mechanical oscillator to non-stationary random excitation. 7. L. L. BUCCIARELLI and C. Ku0 1970 Journal of Applied Mechanics, American Society of Mechanical
Engineers
37, 612-616.
Mean square response of a second order system to
nonstationary random excitation. 8. T. HASSELMAN 1972 Journal of the Engineering Mechanics Division, American Society of Civil Engineers 98, 519-530. Linear response to nonstationary random excitation. 9. R. B. COROTIS and E. H. VANMARKE 1975 Journal of the Engineering Mechanics Division, American Society of Civil Engineers 101,623-637. Time dependent spectral content of system response. 10. J. K. HAMMOND 1968 Journal of Sound and Vibration 7,393-416. On the response of single and multidegree of freedom systems to non-stationary random excitations. 11. D. A. GASPARINI 1979 Journal of the Engineering Mechanics Division, American Society of Response of MDOF systems to non-stationary random excitation. Civil Engineers 105,56-61. 12.D. A. GASPARINI and A. DEBCHAUDHURY (to appear). Dynamic response to non-stationary non-white excitation. The use of 13. M. SAKATA and K. KIMURA 1979 Journal of Sound and Vibration 67,383-393. moment equations for calculating the mean square response of a linear system to non-stationary random excitation. 14. T. K. CAUGHEY 1963 JournaloftheAcousticalSocietyofAmerica 35,1707-1711. Equivalent linearisation techniques. 15. M. SAKATA and K. KIMURA 1980 Journal of Sound and Vibration 73,333-343. Calculation of the non-stationary mean square response of a non-linear system subjected to non-white excitation. 16. S. H. CRANDALL 1963 Journal of the Acoustical Society of America 35, 1700-1705. Perturbation techniques for random vibration of nonlinear systems. 17. J. E. GOLDBERG, J. L. BOGDANOPP and D. R. SHARPE 1964 Bulletin of the Seismological Society of America 54, 263-276. The response of simple nonlinear systems to a random disturbance of the earthquake type. 18. L. D. LUTES and V. S. SHAH 1973 Journal of the Engineering Mechanics Division, American Society of Civil Engineers 99,715-734. Transient random response of bilinear oscillators. 19. Y-K. WEN 1975 Journal of the Engineering Mechanics Division, American Society of Civil Engineers 101,389-401. Approximate method for nonlinear random vibration. 20. R. H. TOLAND and C. Y. YANG 1971 Journal of the Engineering Mechanics Division, American Society of Civil Engineers 97, 791-807. Random walk model for first-passage probability.
28
J. B. ROBERTS 1972 International Journal of Non-Linear random vibration of non-linear structures. J. B. ROBERTS 1977 Journal of Sound and Vibration 50, 145-156. Stationary response of oscillators with non-linear damping to random excitation. W. D. IWAN and P-T. D. SPANOS 1978 Journal of Applied Mechanics, American Society of Response envelope statistics for nonlinear oscillators with Mechanical Engineers 45,170-174. random excitation. P-T. D. SPANOS and W. D. IWAN 1978 Journal of the Engineering Mechanics Division, Computational aspects of random American Society of Civil Engineers 104, 1403-1415. vibration analysis. P-T. D. SPANOS 1976 Report EERL 70-04, Earthquake Engineering Research Laboratory, California Institute of Technology, Pasadena, California. Linearisation techniques for non-
21. R. H. TOLAND, C. Y. YANG and C. K. Hsu Mechanics 7, 395-406.
22. 23. 24. 25.
Non-stationary
linear dynamical systems. 26. P-T. D. SPANOS 1978 International Journal of Solids and Structures 14, 861-867. Nonstationary random vibration of a linear structure. 27. P-T. D. SPANOS 1978 Journal of Structural Mechanics 6, 289-302. Energy analysis of structural vibrations under modulated random excitation.
28. P-T. D. SPANOS 1979 Journal of the Acoustical Society of America 65,404-410.
Hysteretic
structural vibrations under random load. 29. P-T. D. SPANOS 1978 International Journal of Non-LinearMechanics
13,249-259. Stochastic analysis of oscillators with non-linear damping. 30. J. B. ROBERTS 1976 Journal of the Engineering Mechanics Division, American Society of Civil Engineers 102, 851-865.
First passage probability for nonlinear
oscillators.
31. J. B. ROBERTS 1978 in Stochastic Problems in Dynamics, (editor B. L. Clarkson) London: Pitman Publishing Limited. Probability
of first passage failure for lightly damped oscillators. First passage time for oscillators with non-linear restoring forces. J. B. ROBERTS 1978 Journal of Sound and Vibration 60, 177-185. The energy envelope of a randomly excited non-linear oscillator. J. B. ROBERTS 1978 Journal of Applied Mechanics, American Society of Mechanical Engineers 45, 175-180. First passage time for oscillators with non-linear damping. R. L. STRATONOVITCH 1963 Topics in the Theory of Random Noise, Volume 2. New York: Gordon and Breach. N. BOGOLIUBOV and A. MITROPOLSKI 1963 Asymptotic Methods in the Theory of Non-linear Oscillations. New York: Gordon and Breach. S. T. ARIARATNAM 1978 Journal of Applied Mechanics, American Society of Mechanical Engineers 45, 964. Discussion of “Response envelope statistics for nonlinear oscillators with random excitation” by W. D. Iwan and P-T. D. Spanos. R. Z. KHASMINSKII 1966 Theory of Probability and its Applications 11, 390-406. A limit
32. J. B. ROBERTS 1978 Journal of Sound and Vibration 56, 71-86. 33. 34. 35. 36. 37. 38.
theorem for solutions of differential equations with a random right hand side. APPENDIX Consider the conditions under which equation (43) yields a distribution p(v; t) with a which is invariant with respect to time. Let
shape
(Al)
r.4= u/(1 -4),
where q is some function of t only, such that q + 0 as t + co. Then, from equation (43), it is found that the transition density ~(~10; t) is governed by the equation ;
(
u$
*
C(u-qu)
=~I[BB(u-qu)-?rS(wo)l~}+?rs(~~)~ I > (1-q) pI .
W9
If p(uJ0; t) is to be explicitly independent of time it follows that it must also be the stationary density function: i.e., it must satisfy (A3)
TRANSIENT
NON-LINEAR
SYSTEM RESPONSE
29
On subtracting equation (A3) from equation (A2) one has
a ,(~~~)=~{a[B(~-q~)-B(u)$)+rrS(oo)~([~(r~~)-C(~)]~],
(A4)
An inspection of equation (A4) shows that, in general, it cannot be satisfied if q is a function of time only. The only conditions under which this is possible appear to be when both B(v) and C(u) are proportional to L): i.e., B(u) = (YlU,
C(V) = (YZD,
(A5)
where (Y~and (Y~are constants. Equation (A4) then reduces to
which can be satisfied if dqldt = -al/3q.
(A7)
Since p(u ; t) approaches a delta function as t + 0, it is necessary that q + 1 as t + 0 (see equation (Al)). With this initial condition the solution to equation (A7) is simply q = e-+‘* (A8) For the case of the power law oscillator, discussed in section 4.1, B(u) = C(u) = au (i.e., (~2=(w) and equation (A8) agrees with equation (60). The foregoing analysis thus verifies that in this case a time-invariant form of distribution, which is the same as that of the stationary distribution, is possible. This form is given by equation (64). (~1 =
TRANSIENT
RESPONSE
OF NON-LINEAR
TO RANDOM
SYSTEMS
EXCITATION
J. B. ROBERTS School of Engineering
and Applied Sciences, University of Sussex, Falmer, Brighton BN19QT, England (Received 6 June 1980)
An approximate method for computing the transient response of non-linear oscillators to suddenly applied, wide band random excitation is developed. The technique is derived initially for systems with non-linear damping, by using the method of stochastic averaging, and then generalized to deal with oscillators with arbitrary non-linear restoring forces. Predictions from the theory are compared with corresponding digital simulation estimates, in a number of specific cases.
1. INTRODUCTION
In many random vibration problems it is reasonable to model the excitation and response processes as stationary random processes. However, it is important to realize that such modelling is always an approximation, since all real excitation processes are nonstationary to some extent. In some applications there is a very severe degree of nonstationarity in the statistical character of the input which may be best represented as an abrupt change. Examples of such non-stationary excitations are wind gusts, earthquakes, ignition transients in missiles and spacecraft, and aircraft landing loads. In all such cases one is concerned with studying the transient response of structural systems to stochastic excitation. The simplest transient problem, and that which has been most extensively studied, is to compute the statistics of the response of a structural system to suddenly applied, stationary excitation, assuming that the system is initially at rest. For linear systems this problem can be solved fairly easily by a variety of methods, most of which are based on the convolution and Fourier transform relationships (e.g., see references [l-13]). A number of timedependent statistics, such as the mean square response us. time, have been computed by these methods for single degree of freedom systems and for some multi-degree of freedom systems. All real systems are to some extent non-linear, however, and it is thus a matter of practical concern to develop methods of computing the transient response of non-linear systems to stochastic excitation. Unfortunately, it is difficult to extend the standard linear methods to the non-linear case. For the computation of the mean square response Sakata and Kimura have shown that it is possible to extend their method [ 131, for non-stationary excitations, by using the method of equivalent linearization (e.g., see reference [14]), although there are very considerable complications in this approach [ 151. In principle the perturbation method (e.g., see reference [16]) can also be extended to the non-stationary case but again the details are complicated and such an approach has not yet been attempted. Faced with these difficulties, several workers have resorted to the direct digital 11 0022-460X/81/010011
+ 19 $02.00/O
@ 1981 Academic Press Inc. (London) Limited
12
J. B. ROBERTS
simulation approach to obtain estimates of non-linear system response (e.g., see references [17] and [18]). However, when the stochastic excitation is wide band in character an alternative theoretical approach is possible, which is not derived from linear methods. This involves setting up an appropriate Fokker-Planck Kolmogorov (FPK) diffusion equation for the transition density function, and solving this to find the transient response. Unfortunately, even for single degree of freedom systems, the transient response cannot be found analytically and a numerical solution procedure must be used. This approach has been used by Wen [19] to obtain various response statistics for the Duffing oscillator. An alternative numerical procedure has been given by Toland and Yang [20] and Toland et al. [21]. In the FPK approach it is possible to simplify the analysis considerably if the system is a lightly damped oscillator with a small degree of non-linearity. One can then exploit the fact that the amplitude envelope process is approximately a Markov process, governed by a one-dimensional FPK equation (e.g., see reference [22]). This one-dimensional equation is much easier to solve than the exact, two-dimensional FPK equation [19-211. In particular, one can readily find the time-dependent transient response of the amplitude envelope, to a non-stationary wide band input, from the FPK equation. Such results have been obtained for a variety of non-linear oscillators by Iwan and Spanos [23, 241 and Spanos [25-291, using a perturbation technique to find an approximate series solution in terms of eigenvalues and eigenfunctions. It is the purpose of this paper to show that the FPK method, as applied to the amplitude process of a non-linear oscillator, yields considerably more information than just the statistics of the amplitude process itself. Firstly it is demonstrated that it also yields an approximation to the time-dependent joint distribution, p(x, i ; t), of the displacement and velocity processes: from this distribution a wide variety of statistics can be calculated, such as level crossing probabilities. This part of the paper is essentially a generalization of the analysis given in reference [22]. Secondly, it is shown that this approach can be further generalized to deal with oscillators with arbitrary non-linear restoring forces, by using the fact that the energy envelope is Markov [30-331. By solving the FPK equation for the energy envelope it is shown that one can obtain not only the statistics of the energy envelope but also an approximation to the joint distribution p(x, f ; t). For the special case of oscillators with restoring forces governed by a power-law relationship, a simple analytical solution is obtained, which generalizes the previously known result for a linear oscillator. In other cases the appropriate FPK equations are solved by using a random walk analogue [34,32] to form a simple numerical scheme. To assess the accuracy of the approximate theory which is developed here, results obtained from it are compared with corresponding digital simulation estimates, for a variety of cases.
2. THE LINEAR
OSCILLATOR
To introduce the proposed methods it is convenient to begin with the case of a linear oscillator, since here exact solutions are obtainable by standard procedures. The appropriate equation of motion is i+225ooi
+&x
= y(i),
(1)
where n(t) is the displacement response, l is the non-dimensional damping factor and o. is the undamped natural frequency. y(t) will be assumed to be a stationary random process, with zero mean and of wide-band character, which is suddenly applied (i.e., “switched on”) at time t = 0.
TRANSIENT
NON-LINEAR
SYSTEM
13
RESPONSE
The spectrum of y(t) may be defined as
S(o)=& 1-1 W(T)ei‘“’d7,
(2)
where w(7) is the covariance function for y(t): i.e., (3)
w(~)=E{YWY(~+~)I.
It will be assumed that S(o) is smooth in the neighbourhood of w = 00 and that the damping is light. In these circumstances y(t) can be well approximated as a white noise process, n(t), with spectrum S,(w) = S(wo) a constant [l]. The covariance function for n(t) is then W”(7) = 2&(oo)S(7). (4) With y(t) modelled as a white noise it is fairly easy to compute the mean square response of the linear oscillator, CT:,by using the convolution integral and impulse response function (e.g., see reference [l]). The result is
&> -=
1-e-2w
cl2
l+_ [
x2
l-l2
.
2
sm
l@O
_
u’+dl-gz
-sin2fGt
1
(5)
, t>O,
where 6 = 6JoJl
(6)
is the damped natural frequency and CT’is the stationary mean square response, given by Cr2= 7rs(wo)/256J:.
(7)
As 5 -* 0, equation (5) approaches the asymptotic result &(t)/a’
= 1 - e-2yoo’.
(8) In Figures l(a) and (b) r:(t)/ a2 is plotted against the non-dimensional time CN, where N = oot/2r (the number of “cycles”), for [ = 0.05 and O-20, respectively, according to I
Statlanary
(b)
response
Statlonory
I
I
response
0.7 -Jb
0.6
N‘ b’
05 0.4 0.3 0.2 0.1
0
0.1
02
0.3
0
0.1
0.2
O-3
CN Figure 1. Variation of the mean square displacement with time for the linear oscillator. Comparison between the exact solution and the asymptotic solution. (a) [ = O-05; (b) 5 = 0.20.
14
J. B. ROBERTS
equation (5). In both cases a comparison is made with the asymptotic result of equation (8). It is seen that the effect of the sine terms in equation (5) is to produce an oscillation, of frequency 26, superimposed on the smoother variation given by the asymptotic result. The amplitude of this oscillation is of order 5 and the oscillation is seen to be fairly insignificant in the case [ = 0.05. It can be concluded that, for light damping (5 < 0.1, say) of the kind usually encountered in engineering structures, the asymptotic result is a good approximation. This suggests that, in developing a suitable approximate method for lightly damped oscillators, which may not necessarily be linear, one should apply some kind of averaging, or smoothing, procedure, which will eliminate relatively small oscillatory terms such as those appearing in equation (5). Such a “stochastic averaging” procedure has been developed by Stratonovitch [35], and is a generalization of the well-known deterministic averaging method of Bogoliubov and Mitropolski [36]. The averaging method is based on the amplitude envelope process, A(r), where A*(t) = x*(t)+i*(t)/&
(9)
which has a simpler statistical structure than either x(t), or i(t). This can be demonstrated quite easily by considering the mean square amplitude E{A*(t)}. Using the impulse response function again one finds that the mean square velocity, a,*, is given by
af0)
a20Zo=
l-e-2co~r
[
x2
l+-
.
l-l2
2
s1n
iSJO
_
-sin2Gt
“Q-q2
.
I
t>O.
(10)
Hence, from equation (9), w*w1=
I_
1 +
e-2&o,r
2a2
[
x2
sinZ
1_52
&
I *
(11)
The oscillatory term here is of order t*, which will thus be extremely small when the damping is light (cf. equations (5) and (lo), where the dominant oscillatory term is of order 4’). Thus a very good approximation to equation (11) is E{A*(t)} = 2g*(l -e-2rWot)
(12)
and one thus has the asymptotic result (see also reference [25]) a: = &&
= i~{A*(t)}.
(13)
This shows, in a special case, how the smoothed (or averaged) behaviour of the statistics of x(t) and i(t) are related to the statistics of A(t). In the next section a more general relationship between x(t), i(t) and A(t) will be derived, for a class of non-linear oscillators which includes the linear case. 3. THE STOCHASTIC
AVERAGING
METHOD
Consider now a non-linear oscillator, with an equation of motion of the form i+@F(x,k)+&x=y(t),
(14)
where /3 is a coefficient which is small in magnitude. Again y(t) is assumed to be a wide band process, with zero mean. Since the energy dissipation will be small (@small) x(t) and i(t) will have an oscillatory, or narrow band, character and it is meaningful to define an amplitude process, A(t), and a
TRANSIENT
NON-LINEAR
SYSTEM
RESPONSE
15
phase process, 4(t), as follows: x=Acos@,
li=--woAsin@,
(15)
where @=oot+q!L
(16)
A(t), x(t) and i(t)
are then related according to equation (9). The equation of motion can be recast as two simultaneous equations in A(t) and 4(t). On applying the averaging procedure due to Stratonovitch [35] to these equations (see the Appendix in reference [22]) one eventually obtains the approximate equations A = -(p/wo)G(A)
+{~S(wo)/2Ao~}-[nS(wo)]“‘n(t)/wo,
4 = -(PIAwo)H(A)
(17)
- W(oo)11’2N)lAoo,
(18)
F(A cos @, -w. A sin 0) sin ad@,
(19)
F(A cos @, -ooA
(20)
where 2.77
G(A)=-l
2lr I 0
and n(t) and f(t) are mutually independent, such that E{n (t)n 0 + 7)) = WW(t
sin @) cos @ d@,
stationary Gaussian white noise processes,
+ r)] = S(T),
E{n(t)&+T)}=O.
(21722)
It follows from equations (17) and (18) that the joint process A(t), 4(t) is a twodimensional Markov process with a transition density function p(a, c$]u~,c$~;t) governed by the FPK equation
aP a at au
-&-
PWU) ap+dbo) Mwo) %4-y--~P +II au0 a4 -2+7y. Uw. 20~ 0
a2p 1 a2p [ aa
a
a4
3
(23)
(For the special case H(A) = 0, this agrees with the result given in reference [22].) Furthermore, it follows from equation (17) that A(t) is a one-dimensional Markov process, with a transition density function, p(u]uo; t) governed by the FPK equation
ap a -=at
aa
&?(u)-$$$}p]+%$ o. 0
(24)
Here the transient response is sought, with the oscillator assumed to be initially at rest, at time t = 0. The time dependent joint distribution of A(t), 4(t) is then given by
0, 4; t) =p(a, 41% 40; 0,
(25)
where p must satisfy equation (23). On comparing equations (23) and (24) it is found that the axisymmetric form Pb,
4; t)=pb;
tmr,
(26)
where p(u; t) = p(u (0; t) satisfies equation (23). Thus one has the remarkable result that the time dependent solution of equation (24), p(uIO;t), also gives the solution to equation (23). This only applies, of course, with the initial condition A0 = 0.
16
J. B. ROBERTS
From equation (26), by transformation of variables from a, 4 to x, i’, one can express the joint distribution p(x, i; t) of x(t) and i(t) in terms of ~(a; t). The result is P(& i; t) = (claMa;
t),
(27)
where a = (X2+.P/o~)1’2
(28)
and c is a normalization constant. Equation (27) can be used to compute various statistics of x(t) and i(t). For example, the mean square of r(t) is given by -x
dxdi=c
u*p(u ; t) cos* @ dud@ = E{A*(t)}/2,
(2% which is in agreement with equation (13). 3.1. SOLUTION OF THE FPK EQUATION To find ~(a; t), and hence p(x, i; t), the time-dependent required, subject to the initial condition
solution of equation (24) is
p(u; 0) =6(u).
(30)
Unfortunately, a closed form analytical solution seems possible only in the linear case. This well known result is as follows (e.g., see reference [35]): a pb;d=
2
wo(1-4)
exp
(
-
a
2
2&(1-4)
(31)
I
where 4 = exp
cu=a/a
(-25w0tL
W33)
and CTis given by equation (7). Hence, from equation (29), it is readily found that af(t)/a*
= 1 - e-2Loof,
which, as previously shown, is asymptotically correct as [+ 0 (equation equation (31) represents simply a time-dependent Rayleigh distribution, ~(a’; t) = a’ exp (-a’*/2)
(34) (8)). Hence, (35)
where (Y’= a/&&),
(36)
as Spanos has noted [26]. For the non-linear case a numerical solution of equation (24) is necessary. For this purpose the random walk analogue discussed in reference [34] is useful. This involves replacing A(t) by a discrete random walk process, R(ti), which can only assume the amplitudes & = k Su (k = 0, 1, . . .) at time tj = j 8tG = 0, 1, . . .). This process is such that, from state (&-, it can only move to either state (&+I, ti+l) with probability pk, or to state (a&l, ti+l) with probability qk = 1 -pk. With an appropriate choice of pk it can be shown that R(tj) is governed by the same diffusion equation as A(t) (equation (24)), in the limit as St + 0, provided that tj)
@a)* = {7rS(oo)/o~} (at).
(37)
TRANSIENT
NON-LINEAR
17
SYSTEM RESPONSE
In practice the theoretically infinite range of ak must be truncated by choosing some upper limit to k (M say); it is then convenient to set an absorbing barrier at the level UM. The computed distribution, p(& ; ti) will be virtually independent of the value of M, if this is sufficiently large. At ak = 0 the random walk process becomes singular, but this difficulty can be overcome by placing a reflecting barrier at this level. Results obtained from this numerical method will be compared later with corresponding digital simulation results.
4. OSCILLATORS
WITH
ARBITRARY
NON-LINEAR
RESTORING
FORCES
The stochastic averaging method, as described above, is not a useful technique for examining the effect of non-linear restoring forces on the transient response. For example, consider the case of a Duffing oscillator with linear damping; the F(x, i) function in equation (14) is then F(x, i) = i + E&X3.
(38)
On evaluating the G(A) function (see equation (19)) it is found that the cubic term does not contribute, i.e., (39)
G(A) = ooN2,
as for the linear oscillator. Iwan and Spanos have attempted to overcome this difficulty by combining an equivalent linearization method with stochastic averaging, and treating the equivalent frequency as amplitude dependent [23,25,26]. However, they include terms which are not consistent with the order of approximation inherent in the stochastic averaging method and the validity of their approach is doubtful [37]. Here an alternative approach is followed which is based on a consideration of the energy envelope. Consider an oscillator governed by the following equation of motion: f+/?F(x,1)+G(x)=y(t), where G(x) represents an arbitrary process, V(t), is defined as
non-linear
V(r) = P/2 + U(x),
(40)
restoring force. The energy envelope
U(x) = ix G(t) dS.
(4~42)
Jo
Here i2/2 is the kinetic energy of the oscillator and U(x) is its potential energy. In the case where G(x) =&x, previously considered, U(x) = o&*/2 and then V(t) = &A*(t)/2, where A(t) is defined by equation (9). It has been shown previously by the author that V(t) is approximately a one-dimensional Markov process when p is small [30-331. The appropriate FRK equation for the transition density function p(u]uo; t) of V(t) is
ap a -$= ~[BBbd-
d%JoxlP~+
Moo)
gz
(43)
[C(U)Pl,
where 1 B(u) = &A(u)
F(x, J2(v - U)) dx,
Jv-vdx,
(44)
18
J. B. ROBERTS
(45) The integration range, R, in these integrals is such that V(x) < ~1. By solving equation (43), with the initial condition o. = 0, one obtains the timedependent distribution of V(t), p(v; t) =p(vlO; t),
(46)
with the initial form p(u; 0) = S(v). This gives useful information on the rate at which stationarity is achieved. However, one can go considerably further than this by using the fact that, for small p, the joint density of x(t) and V(t), p(x, u; t), is related to p(u; t) according to [30] p(x, u; t) =p(v; t)/2A(o)&
- U.
Upon reverting to x, f variables, equation (47) can be readily transformed expression for the joint density of x(l) and i(t), i.e., p(x, f; t) =p(v;
t)/hA(u).
(47) into an (48)
In the case of linear damping equation (48) agrees with the known exact result for the stationary distribution. Furthermore, in the case where G(x) = wfx then p(u; r) =P(a;
r)l&,
(49)
where tr = w&‘/2, and equation (48) reduces to equation (27). Although equation (27) was derived by using the stochastic averaging procedure, which has a rigorous foundation [38], a correspondingly rigorous derivation of the generalization given by equation (48) seems difficult to achieve. The derivation given in reference [30], whilst physically plausible, is admittedly rather heuristic in nature. A rather obscure argument in support of equation (48) has also been given by Stratonovitch [35]. Some physical insight into the validity of equation (48) can be obtained by considering the phase plane: x, f. In the absence of excitation, and with p = 0 in equation (40), trajectories in the phase plane will be simply orbits corresponding to a fixed value of u =1*/2+ V(x). When the energy dissipation per cycle is small (p small) and the excitation is correspondingly weak, each “cycle” in the phase plane will approximate (in a macroscopic sense) to a constant ~1orbit. For a suddenly applied input the orbit level will gradually increase from zero. In these circumstances it is reasonable to expect that a constant D contour in the phase plane is also a contour of constant probability density, i.e., P(X, i; r) =f(u; 0,
(50)
wheref(v; t) is some function of 0 and t. It follows from equation (50) that the joint density of V(t) and x(t) is given by p(u, x; r) =f(v; t)/_.
(51)
On integrating with respect to x one has p(u;t)=
OJp(v, x; t) dx =f(o; t)&A(v); I -al
(52)
hence (53) Combining equations (50) and (53) one recovers equation (48).
TRANSIENT
NON-LINEAR
SYSTEM
19
RESPONSE
By integration of equation (48) one can find various statistics associated with x(t) and x(t). For example, consider the evaluation of the mth moment of x(t), E{x”(t)}, E{x”(t)}=
Irn Irn x”p(x, x; t) dx dl --g, -ca
x”p(x, u; t) dx dv
=
where
4.1.
ANALYTICAL
SOLUTION
The FF’K equation for V(t) (equation (43)) can be solved analytically in the case where the damping is linear (i.e., F(x, II) = x) and the restoring force has the power-law form
G(x) = k]x]” sgn (x).
(56)
For this purpose it is convenient to define the non-dimensional
energy variable
Z(t) = v(t)/yka”+’
(57)
where (58) and (T is the stationary standard deviation of x(t). Then, as shown in reference [31], the solution to equation (43), when expressed in terms of the transition density p(z]zo; t) of Z(t), is as follows:
l p(zlzo;t>= (l_q)($Jp'(*)exP
{-‘;;_yp’ 51, )Ip(
(59)
where P = (1/a) - 1,
q ; e-olpr
(Y=2(v+l)/(v+3),
(60)
and I,( ) is the modified Bessel function. Here p(z; t) only is required; hence, setting z. in equation (59) to zero and using the fact that I,(x)~{l/~(l+p)}(x/2)p
asx+O
(61)
it is found that p(z; t> =
1 r(1 +p)(l -q)‘+‘zp
-_ f
exp (
l-q > -
(62)
In the linear case p = 0, (Y= y = 1, /3 = 24’~~ and equation (62) is equivalent to equation (31). It is interesting to note that the shape of the distribution of Z(t), as given by equation (62), is invariant with respect to time. Thus, if z’ = z/(1 -4)
(4
then equation (62) can be written as p(z'; t)= {i/r(i+p)}(~')p
exp
(-z').
(64)
20
J.B. ROBERTS
In the special case where t + co, then q + 0, z'= z,and equation (64) gives the stationary density function. This conclusion is a generalization of that given earlier for the particular case of a linear oscillator (see equations (31), (35) and (36)). The invariance with respect to time of the form of the distribution of energy, found here for the case of the power law oscillator, is due to the fact that B(n) and C(U), in equation (43), are identical and proportional to U. It is shown in the Appendix that time-invariance of the shape of p(v; t) occurs only when B(u) and C(v) are both proportional to ZI. The moments of Z(t) can be easily computed from the above results by suitable integration. It is found that E{Z”(t)}=(l-q)(p+m)(p+m-1).
. . (p+l),
(65)
and in the special case where m = 1 this reduces simply to E{Z(t)) = (1 -q)(p + 1) =
(1--4)/a.
(66)
In Figure 2(a) the non-dimensional mean energy E{V(t)}/k~‘+” = yE{Z(t)} is plotted against the non-dimensional time, pt, for Y = 0, 1, 2 and 3, according to equation (66).
Figure 2. (a) Variation of the mean energy with time for the power-law oscillator; effect of power index; (b) variation of the mean square displacement with time for the power-law oscillator; effect of power index.
Similarly the moments of x(t) can be found by combining equations (54) and (62), and integrating. For example, the mean square response of x(t) is found to be given by the following remarkably simple formula: &(t) = a*(1 - q)2'(1+v).
(67)
For the linear case (Y = 1) this result reduces to equation (34). Figure 2(b) shows the variation of &(t)/a* with fit for y = 0, 1, 2 and 3, according to equation (67). 4.2. NUMERICAL SOLUTION OF THE FPK EQUATION In general equation (43) cannot be solved analytically and a numerical solution procedure is required. For the results to be presented in the next section an extension of
TRANSIENT
NON-LINEAR
SYSTEM
RESPONSE
21
the random walk analogue has been used, as described in reference [32]. Briefly, the method is based on a non-linear co-ordinate transformation, which transforms equation (43) into a diffusion equation of the same form as a limiting random walk (at + 0). The co-ordinate transformation must be performed numerically but need be carried out only once, prior to the numerical diffusion process. The time step and amplitude increment must be related, in a similar way to that indicated in equation (37).
5. COMPARJSON
WITH
SIMULATION
RESULTS
To test the approximate theory developed a number of comparisons have been made between theoretical predictions and simulation estimates, as follows. 5.1,
SIMULATION
METHOD
The simulation method consists of generating a number of realizations of the transient response, from t = 0 to t = T, by digital means. For this purpose the input process y(t) is approximated as a step function which jumps from one level to another at equi-spaced time instants, at spacing dt. The levels are proportional to a series of independent, Gaussian random numbers. Such a step process is a very close approximation to a white noise process when At is small. The equation of motion is integrated by using a standard fourth order Runge-Kutta algorithm, with an integration step length of At. To obtain reasonably stable estimates of parameters such as the mean energy and mean square displacement, 4000 realizations were generated, for every case. These realizations were ensemble averaged to provide the required statistical estimates. 5.2.
OSCILLATORS
WITH
NON-LINEAR
DAMPING
It was shown earlier that, for the case of linear damping, the stochastic averaging theory gives a result for the mean energy envelope which has an error of order [*, whereas for the mean square displacement the error is of order f: One can thus anticipate that, in the more general, non-linear case the approximate theory is likely to give more accurate estimates of total energy statistics than of displacement (or velocity) statistics. To test the theory for non-linear damping an oscillator with a cubic damping is considered-the equation of motion is as follows: ~+2~~o(~+E1~3)+W~X=y(t).
(68)
If u. is the standard deviation of the stationary response in the linear case (&r = 0) then equation (68) can be conveniently non-dimensionalized in terms of X=x/u,
and
r =oot.
(69)
It then becomes X+25(X+eEX3)+X=
Y(r),
(70)
Y(r) = Ywbo&.
(71)
where &T= “;I_&, The corresponding
non-dimensional
energy is then v* = x2/2 +X2/2.
(72)
In Figure 3(a) the theoretical variation of E{ V*} with non-dimensional time cN(N = oot/27r) is shown, for ET =0, O-2 and 1.0. These result? were obtained by using the random walk method to solve equation (24) numerically. Also shown in this figure are
22
J. B. ROBERTS
0
1
I
I
/
O-04
0.06
042
M6
I
SN
Stationary ---------------
response
(
0.6 05 0.4 0.3 0.2 0.1
0
I
-+-
I
4
N Figure 3. (a) Variation of the mean energy with time for an oscillator with a cubic damping term, comparison of theory with simulation results; simulation estimates: [ = 0.04: v, ET = 0; X, ET = 0.2; +, ET = 1.0; (= 0.10: 0, ET = 0; 0, B: = 0.2; 0, ET = 1.0; (b) variation of the mean square displacement with time for an oscillator with a cubic damping term, C = 0.04, comparison of theory with simulation results; theory, -, simulation; (c) variation of the mean square displacement with time for an oscillator with a cubic damping term, 5 = 0.10, comparison of theory with simulation results; -, theory, m, simulation.
corresponding simulation estimates of E{ V*}, for 4’= 0.04 and 4’= O-10. The time step, AN = woAt/27r, in the simulation was 0.02. There is good agreement between the theoretical curves and the simulation estimates at all three values of ET. Figure 3(b) shows the computed variation of the mean square of X(r), (T;, with N for 4’= 0.04. It is noted that these theoretical curves can be obtained directly from Figure 3(a) since, according to the stochastic averaging theory, r&(f) = E{ V*} (see equation (29)). The corresponding simulation results vary in an oscillatory way with time, as one might expect in view of the exact linear result (see Figures l(a) and (b)): to bring out this characteristic the
TRANSIENT
NON-LINEAR
SYSTEM
RESPONSE
23
simulation results are here plotted as a continuous line, rather than as a series of points, as in Figure 3(a). Again it is noted that good agreement is obtained between the two sets of results. One would expect the exact solution for ET > 0 to be oscillatory, as in the linear case, and this is indicated by the simulation results, although this behaviour is to some extent obscured by the inherent statistical uncertainty. A similar comparison is made in Figure 3(c) for the case 5 = 0.10. There is a tendency for the amplitude of the oscillations in the simulation results to be greater at the higher damping factor: again one can anticipate this by reference to the linear case and also from the nature of the stochastic averaging theory. 5.3.
OSCILLATORS
WITH
NON-LINEAR
STIFFNESS
5.3.1. Power-law oscillators with linear damping Again it is convenient at the outset to non-dimensionalize is
the equation of motion, which
i+pi+k(xj”sgn(x)=y(t).
(73)
If X=x/u
and
(74)
r’=Pt
then equation (73) can be written as X+X If a non-dimensional
+ (k~“-‘/~2)~X~” sgn (X) = y(t)/&.
(75)
damping coefficient is defined as x = /3/2[ka”-‘1”’
(76)
(in the linear case x = l), then it is evident that the non-dimensional response depends only on Y and x. Figure 4(a) shows the variation of the mean non-dimensional energy yE{Z(t)) with time for Y = 0 and v = 1 (see also Figure 2(a)). These were computed from equation (66). Also shown are simulation results for ,y = 0*04,0*10 and O-50. At the highest of these damping values there is some discrepancy between the simulation results and the approximate theoretical results, but the agreement is much better at the lower values of damping. Figures 4(b) and (c) show comparisons between theory and simulation estimates for the mean square displacement, a& us. time, for x = 0.04 and x = O-10, respectively. It is clear that the agreement is best at the lower value of damping as one would expect. The oscillations in the simulation results are only partly due to statistical variability: there is little doubt that the exact solution would be oscillatory, the ripple having a frequency and amplitude similar to that exhibited by the simulation results. 5.3.2. Oscillators of the Dufing type Finally an oscillator of the Duffing type has been considered; i.e., one with an equation of motion as follows: i +2&k Again, non-dimensionalization one obtains
+0&L +&2x3) = y(t).
in terms of X and T, as given by equations (69) is useful: X+21X+X+&X3=
where Y(7) = y(t)/ao&
(77)
ET = &$EZand differentiation
Y(t),
is now with respect to T.
(78)
24
J. B. ROBERTS
I.0
0.9
0.8 0.7 0.6
N
2
0.5 0.4 0.3 0.2 01
0
0.1
0.2
03
0
0.2
03
04
P, z
Figure 4. (a) Variation of the mean energy with time for the power-law oscillator, comparison of theory with simulation results; simulation estimates: Y = 0: ?? , x = 0.04; +, x = 0.10; A, * = 0.50; Y = 3:0, x = 0.04; X, x = O-10; V, x = 0.50; - - -, stationary response; (b) variation of the mean square displacement with time for the power-law oscillator ,y = 0.04, comparison of theory with simulation results; (c) variation of the mean square displacement with time for the power-law oscillator, x = 0.10, comparison of theory with simulation results.
In Figure 5(a) the variation of the mean of the total non-dimensional v* = R2/2 +x2/2
+ E:x4/4,
energy, (79)
with time is plotted for E; = 0,0*2 and 1.0. These were computed from the FPK equation for V(t) (equation (43)), by using the extended random walk analogue method. Simulation estimates for E{ V*}, for 5 = 0.04 and O-10, are seen to be in good agreement with the theoretical curves. It is interesting to note that the total energy is only weakly dependent on the extent of the non-linearity. Figures 5(b) and (c) show comparisons between theory and simulation estimates for the mean square displacement, aat), for [= 0.04 and [ = 0.10: respectively. Again the
TRANSIENT
(b) I.0
NON-LINEAR
I
I
I
I
I
2
3
4
25
SYSTEM RESPONSE
_-___-__---___--__----------
0.9
0.8 D7 0.6 “2
0.5 0.4 0.3 0.2 @I
, 0
I I
I
I
I
2
3
4
N
Figure 5. (a) Variation of the mean energy with time for the Duffing oscillator, comparison of theory with simulation results; simulation estimates: 5 = 0.04: 0, ef = 0.2; 0, ei: = 1.0; 5 = O*lO: X, e: = 0.2; +, E; = 1.0; (b) variation of the mean square displacement with time for the Dulling oscillator, C = 0.04, comparison of theory with simulation results; - - -, stationary response; -, theory; MA,simulation; (c) variation of the mean square displacement with time for the Duffing oscillator, c = 0.10, comparison of theory with simulation results. - - -, stationary response; -, theory; m, simulation.
simulation
results
indicate
an oscillatory
solution
which
is close to the present
theory
when
the damping is light. 5.4. ACCURACY OF THE COMPUTED RESULTS It is interesting to enquire briefly into the accuracy of the computed theoretical results, obtained by using the random walk analogue. For this purpose it is convenient to examine some results obtained for the Duffing oscillator; for this case the stationary mean square
26
J. B. ROBERTS
AX Figure 6. Variation of the percentage error in the mean square stationary displacement, as computed by the random walk analogue, with the amplitude increment.
displacement can be computed exactly and compared with the asymptotic result obtained from the numerical diffusion process, as the elapsed time becomes large. In the extended random walk method the range for X(t) is,.initially, divided up into a number of equi-spaced amplitude increments, dx say. Thus dx is a measure of the coarseness of the discretization in the random walk method. Figure 6 shows the variation of the percentage error in the stationary solution, with dx, for Ez = 0.20 and E: = 1.0. This shows that, in both cases, there is an almost linear variation of error with Ax. In the results given earlier a value of Ax = O-01 was chosen; it is clear that the error is then only a few percent, at most. As a final point, it is noted that one can use the known, exact result for the stationary response to correct the results obtained by the numerical diffusion procedure. Of course, I
06
I
I
I
I
-
_.-.------. 0.5 -
04N
c
0
02
c-4
06
k0
Figure 7. Variation of the mean square displacement with time for the Dutling oscillator. ET = 1.0. Effect of correcting the results from the random walk method by using the known exact result for the stationary case. corrected estimate, Ax = 0.01. -* -, Original estimate, Ax = O-07; - - -, corrected estimate Ax = 0.07; -,
TRANSIENT NON-LINEAR SYSTEM RESPONSE
27
for the transient response such a correction will be an approximation. Figure 7 shows the computed variation of crgwith [N for ~5 = 1.0 and an amplitude increment Ax = O-07. For the stationary response the error is 19.6% (see Figure 6). If the ordinates are reduced everywhere by 19.6% one obtains a curve which is in very close agreement with the corrected result for Ax = O-01 (error 2.8%).
REFERENCES 1. T. K. CAUGHEY and H. J. STUMPF 1961 Journal ofApplied Mechanics, American Society of Mechanical Engineers 28, 563-566. Transient response of a dynamic system under random excitation. 2. R. L. BARNOSKI and J. R. MAURER 1969 Journal of Applied Mechanics, American Society of Mechanical Engineers 36, 221-227. Mean-square response of simple mechanical systems to nonstationary random excitation. 3. R. L. BARNOSKI and J. R. MAURER 1973 Journal of Applied Mechanics, American Society of Mechanical Engineers 40, 73-77. Transient characteristics of simple systems to modulated random noise. 4. J. B. ROBERTS 1971 Journal of Sound and Vibration 14,385-400. The covariance response of linear systems to non-stationary random excitation. The covariance response of 5. J. B. ROBERTS 1971 Journal of Sound and Vibration 17,299-307. linear systems to locally stationary random excitation. p. T. K. HASSELMAN 1970 Journal of Applied Mechanics, American Society of Mechanical 189. A comparison of solutions for the mean square response of a simple Engineers 37,1187-l mechanical oscillator to non-stationary random excitation. 7. L. L. BUCCIARELLI and C. Ku0 1970 Journal of Applied Mechanics, American Society of Mechanical
Engineers
37, 612-616.
Mean square response of a second order system to
nonstationary random excitation. 8. T. HASSELMAN 1972 Journal of the Engineering Mechanics Division, American Society of Civil Engineers 98, 519-530. Linear response to nonstationary random excitation. 9. R. B. COROTIS and E. H. VANMARKE 1975 Journal of the Engineering Mechanics Division, American Society of Civil Engineers 101,623-637. Time dependent spectral content of system response. 10. J. K. HAMMOND 1968 Journal of Sound and Vibration 7,393-416. On the response of single and multidegree of freedom systems to non-stationary random excitations. 11. D. A. GASPARINI 1979 Journal of the Engineering Mechanics Division, American Society of Response of MDOF systems to non-stationary random excitation. Civil Engineers 105,56-61. 12.D. A. GASPARINI and A. DEBCHAUDHURY (to appear). Dynamic response to non-stationary non-white excitation. The use of 13. M. SAKATA and K. KIMURA 1979 Journal of Sound and Vibration 67,383-393. moment equations for calculating the mean square response of a linear system to non-stationary random excitation. 14. T. K. CAUGHEY 1963 JournaloftheAcousticalSocietyofAmerica 35,1707-1711. Equivalent linearisation techniques. 15. M. SAKATA and K. KIMURA 1980 Journal of Sound and Vibration 73,333-343. Calculation of the non-stationary mean square response of a non-linear system subjected to non-white excitation. 16. S. H. CRANDALL 1963 Journal of the Acoustical Society of America 35, 1700-1705. Perturbation techniques for random vibration of nonlinear systems. 17. J. E. GOLDBERG, J. L. BOGDANOPP and D. R. SHARPE 1964 Bulletin of the Seismological Society of America 54, 263-276. The response of simple nonlinear systems to a random disturbance of the earthquake type. 18. L. D. LUTES and V. S. SHAH 1973 Journal of the Engineering Mechanics Division, American Society of Civil Engineers 99,715-734. Transient random response of bilinear oscillators. 19. Y-K. WEN 1975 Journal of the Engineering Mechanics Division, American Society of Civil Engineers 101,389-401. Approximate method for nonlinear random vibration. 20. R. H. TOLAND and C. Y. YANG 1971 Journal of the Engineering Mechanics Division, American Society of Civil Engineers 97, 791-807. Random walk model for first-passage probability.
28
J. B. ROBERTS 1972 International Journal of Non-Linear random vibration of non-linear structures. J. B. ROBERTS 1977 Journal of Sound and Vibration 50, 145-156. Stationary response of oscillators with non-linear damping to random excitation. W. D. IWAN and P-T. D. SPANOS 1978 Journal of Applied Mechanics, American Society of Response envelope statistics for nonlinear oscillators with Mechanical Engineers 45,170-174. random excitation. P-T. D. SPANOS and W. D. IWAN 1978 Journal of the Engineering Mechanics Division, Computational aspects of random American Society of Civil Engineers 104, 1403-1415. vibration analysis. P-T. D. SPANOS 1976 Report EERL 70-04, Earthquake Engineering Research Laboratory, California Institute of Technology, Pasadena, California. Linearisation techniques for non-
21. R. H. TOLAND, C. Y. YANG and C. K. Hsu Mechanics 7, 395-406.
22. 23. 24. 25.
Non-stationary
linear dynamical systems. 26. P-T. D. SPANOS 1978 International Journal of Solids and Structures 14, 861-867. Nonstationary random vibration of a linear structure. 27. P-T. D. SPANOS 1978 Journal of Structural Mechanics 6, 289-302. Energy analysis of structural vibrations under modulated random excitation.
28. P-T. D. SPANOS 1979 Journal of the Acoustical Society of America 65,404-410.
Hysteretic
structural vibrations under random load. 29. P-T. D. SPANOS 1978 International Journal of Non-LinearMechanics
13,249-259. Stochastic analysis of oscillators with non-linear damping. 30. J. B. ROBERTS 1976 Journal of the Engineering Mechanics Division, American Society of Civil Engineers 102, 851-865.
First passage probability for nonlinear
oscillators.
31. J. B. ROBERTS 1978 in Stochastic Problems in Dynamics, (editor B. L. Clarkson) London: Pitman Publishing Limited. Probability
of first passage failure for lightly damped oscillators. First passage time for oscillators with non-linear restoring forces. J. B. ROBERTS 1978 Journal of Sound and Vibration 60, 177-185. The energy envelope of a randomly excited non-linear oscillator. J. B. ROBERTS 1978 Journal of Applied Mechanics, American Society of Mechanical Engineers 45, 175-180. First passage time for oscillators with non-linear damping. R. L. STRATONOVITCH 1963 Topics in the Theory of Random Noise, Volume 2. New York: Gordon and Breach. N. BOGOLIUBOV and A. MITROPOLSKI 1963 Asymptotic Methods in the Theory of Non-linear Oscillations. New York: Gordon and Breach. S. T. ARIARATNAM 1978 Journal of Applied Mechanics, American Society of Mechanical Engineers 45, 964. Discussion of “Response envelope statistics for nonlinear oscillators with random excitation” by W. D. Iwan and P-T. D. Spanos. R. Z. KHASMINSKII 1966 Theory of Probability and its Applications 11, 390-406. A limit
32. J. B. ROBERTS 1978 Journal of Sound and Vibration 56, 71-86. 33. 34. 35. 36. 37. 38.
theorem for solutions of differential equations with a random right hand side. APPENDIX Consider the conditions under which equation (43) yields a distribution p(v; t) with a which is invariant with respect to time. Let
shape
(Al)
r.4= u/(1 -4),
where q is some function of t only, such that q + 0 as t + co. Then, from equation (43), it is found that the transition density ~(~10; t) is governed by the equation ;
(
u$
*
C(u-qu)
=~I[BB(u-qu)-?rS(wo)l~}+?rs(~~)~ I > (1-q) pI .
W9
If p(uJ0; t) is to be explicitly independent of time it follows that it must also be the stationary density function: i.e., it must satisfy (A3)
TRANSIENT
NON-LINEAR
SYSTEM RESPONSE
29
On subtracting equation (A3) from equation (A2) one has
a ,(~~~)=~{a[B(~-q~)-B(u)$)+rrS(oo)~([~(r~~)-C(~)]~],
(A4)
An inspection of equation (A4) shows that, in general, it cannot be satisfied if q is a function of time only. The only conditions under which this is possible appear to be when both B(v) and C(u) are proportional to L): i.e., B(u) = (YlU,
C(V) = (YZD,
(A5)
where (Y~and (Y~are constants. Equation (A4) then reduces to
which can be satisfied if dqldt = -al/3q.
(A7)
Since p(u ; t) approaches a delta function as t + 0, it is necessary that q + 1 as t + 0 (see equation (Al)). With this initial condition the solution to equation (A7) is simply q = e-+‘* (A8) For the case of the power law oscillator, discussed in section 4.1, B(u) = C(u) = au (i.e., (~2=(w) and equation (A8) agrees with equation (60). The foregoing analysis thus verifies that in this case a time-invariant form of distribution, which is the same as that of the stationary distribution, is possible. This form is given by equation (64). (~1 =