Reliability of hysteretic systems subjected to white noise excitation

Structural Safety, 8 (1990) 369-379

369

Elsevier

RELIABILITY OF HYSTERETIC S Y S T E M S SUBJECTED TO WHITE NOISE EXCITATION * S.R.K. Nielsen, K.J. Merk and P. Thoft-Christensen Institute of Building Technology and Structural Engineering, University of Aalborg, Sohngaardsholmsvej 57, DK- 9000 Aalborg, Denmark

(Accepted December 1989)

Key words: stochastic differential equations; hysteretic systems; damage indicators; filtered white noise excitation; non-linear polynomial expansion; cumulant-neglect closure; modified cumulant-neglect closure.

ABSTRACT A single-degree-of-freedom bilinear hysteretic system driven by filtered white noise is analysed by means of stochastic differential equations. Failure of the system takes place by first passages at critical levels of so-called damage indicators, which are non-decreasing response quantities measuring at a macroscopic scale the strength and stiffness deterioration. For illustration the cumulative plastic deformation ratio has been used as damage indicator. The white noise excitation is filtered through a time-invariant rational second-order filter. The analytical technique applied is based on a replacement of the original system with non-linear and non-analytical right hand sides for the constitutive equation and damage indicator differential equation by an equivalent system, where these equations are given by a cubic polynomial expansion. The expectations appearing in the equations for the coefficients of the equivalent polynomials are calculated by a multi-dimensional Gram-Charlier expansion. In order to close the infinite hierarchy of moment equations the conventional cumulantneglect closure scheme and a modification with due consideration to the finite probability @yielding are considered. The analytical results obtained are compared to those obtained by Monte Carlo simulation, from which it is concluded that substantial improvements are obtained by applying the modified closure scheme.

* Paper presented at the Euromech 250 Colloquium on Nonlinear Structural Systems under Random Conditions, Como, Italy, June 19-23, 1989. 0167-4730/90/$03.50 © 1990 - Elsevier Science Publishers B.V.

370

1. INTRODUCTION Hysteretic oscillators are mechanical systems in which the dynamic restoring forces are non-conservative with a complicated dependence on the deformation history. For this reason the constitutive equations for the restoring forces must be incrementally formulated. Due to the fact that the incremental stiffnesses are different for loading and unloading, the right hand sides of the constitutive differential equations become non-analytical. During hysteretic deformations the micro-structure of the material is partly damaged due to dislocation migration, development of micro-cracks, etc. In a structural model these damages must be related to the state variables of the system. Differential models for a number of such damage indicators have been formulated for single-degree-of-freedom (SDOF) systems [1,2], as well as for multi-degree-of-freedom (MDOF) systems [3]. From the physics of the problem it is clear that any damage indicator will be a non-decreasing function with time. As is the case for the constitutive equations the fight hand sides of the differential equations specifying the development of the damage indicators may also be a non-linear and non-analytical function of the state variables. Stochastic response of hysteretic structures can be studied based on stochastic differential equations [2-4]. When the loading process is modelled as a white noise process, filtered through a time-invariant rational filter, the state vector made up of state variables of the structural system, the filter and damage indicators becomes a Markov diffusion vector process. The time development of the state vectors is specified by a system of It6-differential equations. The associated F o k k e r - P l a n c k - K o l m o g o r o v equation can hardly be integrated, neither analytically nor numerically, except in the simplest cases. Caughey [5] developed an approximate method of equivalent non-linear differential equations in order to obtain the steady-state response of a class of oscillators with non-linear damping subjected to stochastic excitations. For instance, for an ideal elastic-plastic oscillator subjected to stationary white noise excitation, the variance response is of the non-stationary Brownian motion type. For such hysteretic systems this method is obviously not applicable. Because of these difficulties one has concentrated on obtaining the time-dependent zero-timelag statistical moments of the joint probability density function. Differential equations for these quantities are easily derived, applying the It6-formula to products of the state variables of ascending order, and performing the expectation of the result. The expectations on the right hand side of these differential equations can only be calculated by means of an approximate joint probability density function (pdf), specified with due consideration to the physics of the problem, and with a number of parameters calibrated to the statistical moments available. In cases of even a moderate number of state variables this may lead to serious numerical problems. Equivalent linearization remains the most efficient method, at least for problems of high dimensionality. The idea of the method is to replace the drift vector of the original system by a sequentially updated linear drift vector. The first formulations of the method assumed zero-mean response. Later estimates to non-zero-mean problems were due to Spanos [6] and Baber [7]. Equivalent linearization fails to provide reliable estimates of the variance response for highly non-linear problems. As an example the method predicts a limiting stationary state for an elasto-plastic subjected to stationary white noise excitation in contradiction to the non-stationary Brownian motion type of variance drift observed by numerical simulations [3,8]. Instead it can be attempted to replace the original system with non-linear and non-analytical drift vector by a sequentially updated equivalent system, in which the drift vector is given by an

371

M t h order polynomial expansion of the state variables. If the coefficients of the equivalent polynomium are determined by a least mean square procedure it can be shown that the original and the equivalent systems determine identical mean values and joint moments up to and including order M + 1 provided the same approximate pdf is applied to both systems for the evaluation of expectations entering both systems [3,8]. The polynomial expansion technique was formulated by Nielsen et al. [8] for the zero-mean problem, and by Mork [3] for the non-zero-mean problem. The infinite hierarchy of moment differential equations related to these methods was closed by a cumulant-neglect closure scheme. The objectives of the present study are to present some further developments of the analytical techniques developed by the authors. Especially a modification of the cumulant-neglect closure scheme is suggested, which takes into consideration the finite probability of yielding, formally appearing as delta functions in the joint pdf of the state variables. The analytical results obtained are compared to those obtained by Monte Carlo simulation, and it is concluded that a significant improvement is obtained compared to the conventional cumulant-neglect closure scheme. 2. STOCHASTIC DIFFERENTIAL EQUATIONS FOR SDOF HYSTERETIC OSCILLATOR UNDER FILTERED WHITE NOISE EXCITATION The equation of motion of an SDOF hysteretic oscillator can be written 5~+ 2~0~ + t o 2 ( ( 1 - a ) x + e t z ) = a a ~ + a2 u

(1)

where x and :t are the displacement and velocity of the system, respectively, and z is the hysteretic component of the restoring force, x and z have been made non-dimensional with respect to the elastic limit displacement on the skeleton curve, e.g. x = 1 corresponds to yielding (cf. Fig. 1). o~ is the cyclic linear undamped eigenfrequency, ~ is the damping ratio and a is the post- to preyielding stiffness ratio.

z

a)

~ : ( ~ + (1 - a)z)

b)

X

dx

X B

~ k=0 \

k=0 __ /

/

k=0

~lx'

E k=i/ o/ \k=l

~X cCO

\k=0

} (02otdx

co2

k=l

2 1

TM

D Fig. 1. Bilinear restoring force as a function of non-dimensional displacement x. (1) Non-dimensional hysteretic component. (b) Total restoring force ~2(ax + ( 1 - a)z).

372 The loading of the system is given by aafi + a2u , where a I and a 2 are constants, and u is obtained as the response of a linear filter subjected to a Gaussian white noise excitation. The state vector made up of the state variables of the oscillator, filter and damage indicator is specified by the following system of first-order It6-differential equations

dX(t)=d(X)dt+e(t)dW(t), d ( X ) = A ( X tX

x(t)=

x 2 z xp , u

X(to)=X o

(2) (2a)

0 0 0

e(t) =

0

(2b)

0

a

B(t)

0

1

--aW 2

--2~'W

0

A(X)=

t>to,

--(1

z)

o

0 0

0 0

0

0

0

0

-- a ) 6 9 2

0

a 2

a1

0

0

0

0

o

o

o

o

0 0

0 0

0 _~2

1 -2rl~2

(2c)

X 0 is the initial conditions of the stochastic vector X at the time t = t 0. Zero start will be assumed, i.e. X 0 = 0. The third equation in (2) represents the constitutive equation for the hysteretic component z in differential form. k(2, z) may be interpreted as a state-dependent stiffness. In the present study a bilinear oscillator is considered in which case k(2, z) is given by k(2, z)= 1 -[1 - H(-2)]H(z-

1 ) - [1 - H ( 2 ) ] H ( - z -

1)

(3)

where H ( . ) is the Heaviside unit step function defined by H(x)=

{1, O,

x>~0 x<0

(41

The value of k(2, z) on various branches of the x - z curve has been marked in Fig. l(a). Xp is the non-dimensional accumulated numerical plastic deformation which makes up the damage indicator considered. The differential equation for Xp, as specified by the fourth equation in (2), becomes [1,2]

2p=(1-a)[H(2)H(z-

1)-H(-2)H(-z-

1 ) ] 2 = h ( 2 , z)2

(5)

Equation (5) is easily derived from the sketch shown in Fig. l(b). It follows that the functions k(2, z) and h(2, z) fulfill the symmetry conditions

k(2, z ) = k ( - 2 , - z ) h(2, z ) = - h ( - 2 , -z)

(6) (7)

Provided X 0 = 0, the condition (6) implies E[2] -- E[z] = 0 and Elk(2, z)2] = 0. The fifth and sixth equations in (2) specify the response of the filter. ~2 and ~/ are the cyclic eigenfrequency and damping ratio of the filter, respectively. The considered filter is time-invariant, linear and of second order. Time-varying, non-linear or higher-order filters may be implemented without any difficulties.

373

p(t) specifies the intensity of the white noise. This is assumed in the form

TO is the time of maximum intensity &, and b is a limiting decay rate of the excitation. &, is specified so that the stationary standard deviation of the response of the corresponding linear system under stationary excitation becomes equal to y. &, is then given by

PO=

2b,( b32+ b$, - b,b,b,)

y

(9)

(b, - brb,)

i

b, = 21$ + 23w,

b, = Q2 + w2 + 437jwQ

b, = 23wQ2 + 022q!&

b4 = w2Q2

Finally, { W(t), t E [t,, co]} signifies a unit Wiener process. This is a Gaussian process with the incremental properties E[dW( t)] = 0,

E[dJJ’(t,)

dWt,)l

= (“,;

1

3. DIFFERENTIAL EQUIVALENT

00)

;’ : :’

9

EQUATIONS FOR JOINT CUBIC SYSTEMS

2

CENTRAL

MOMENTS

OF ORIGINAL

AND

Let g( X, t) be an arbitrary function with continuous partial derivations with respect to x and Making use of the It6 differential rule for Markov diffusion processes, and taking the expectation, the following differential equation for the expectation E[ g( X, t)] can be derived [9] t.

-&E[g(X,

t)]

=E[kg(X,

t)]+E[di(X)&g(Xv

t~]++eiejE[&g(Xv I

t)] I

J

(11)

di and e, signify components

of the drift vector d and diffusion term e in eqn. (2). In eqn. (11) and in the following the summation convention has been used over dummy indices. Let

g(x,

t) =

n=l

xi?

qx2... xi”., n> 1 ( where [X,, X2, X,, X,, X5, X6] = [x, i, z, xp, U, zi] and

(12)

x=xi-pi

(13)

Pi=E[Xi] Then

(14)

fii = E[ di]

(-4

fiii = E[ d:XJ] + E[ d;x,“] + eiej

(15b)

j..iijk= E[ dfX;z]

+ E[ d;q$]

/.iijkl = E[ d;XJ~~] +

eiejPkr

+

+ E[ d;Xi”q]

+ E[ d;qzz] eiekPjf

+

eie#jk

054

+ E[ d;xqx] +

ejekPu

+

ejePik

+ E[ d;xXJx] +

eke#ij

(154

374 where d~(X) = d~(X) - E[ d,(X)] /~,,..

.i,,=E[X,, c

.

.

.

(16)

X i,,] c

(17)

Instead of the original system an equivalent non-linear system is considered for which the centralized drift vector is given in terms of a cubic expansion in the centralized state variables [3]

di,eq(X) c

=

Ai + BimS m + CimnX~nS~c c

c

(18)

c c

The idea is, that the joint moments of the two systems are alike in case of a proper choice of the expansion coefficients. In this relation the following proposition can be proved: In case an equivalent polynomium of order M is applied, the propagation of the mean value and joint central moments up to and including order M + 1 will be identical for the original and the equivalent system, provided the same approximate pdf is applied to both systems for evaluation of the expectations including unprovided joint central moments of order M + 1.... , 2 M appearing in the equations for the equivalent coefficients, and provided the coefficients of the polynomial expansion of the equivalent system are determined from a least mean square criterion [3,8]. Obviously this proposition stresses the importance of using the least mean square criterion for the determination of the equivalent polynomial expansion coefficients. The differential equations for the joint central moments of the equivalent system are obtained by inserting (18) into (15b)-(15d) ~tij ~-- Bim~m j + Bjm~m i + Cirnnbtmnj "+ Cjmn~mn i + Dimnpl~mnpj + Ojmnp~mnp i + eiej

(19a)

12ijk = Ail.Lj k ._}_Ajldik + Ak~i j + Biml~mjk + Bjm~mi k qt_ Bkml.tmi j + Cimn#mnjk + Cjmn~mnik + Ckmn~mnij + DimnpP, mnpj k + Djmnp~mnpik + Dkmnp~mnpij ~ijkl = AilXjkl + Ajl'~ikl + Akt-tijl + Alltijk + niml'tmjkl + njm~mikl -']- Bkml'Lmijl +

(19b)

B~,,,l~,,,ij~,

+ Cimn~mnjk l + Cjmn~mnik l'~ Ckrnn~mnijl + Clmn~mnij k + OimnplAmnpjkl + Ojmnp~rnnpikl "t- Dkmnpt.tmnpijl + Dlmnp~mnpij k + eiejlz m + eiekFjt + eied~jk + ejekl~il + ejel~ik + ekedxij

(19c)

If the requirements of the proposition are fulfilled, eqn. (19) results in the same joint central moments as in eqn. (15). The main problem is, that the joint central moments/*,,,ijk and #m,,p,jk of the 5th and 6th order are not provided by the system. Deviation at the solutions occurs if these unprovided moments are calculated by some approximate scheme. In this paper the error in these approximations is investigated.

4. CUBIC EXPANSION TECHNIQUE

In (2) the right hand sides of the third and fourth equations are the only non-linear and non-analytical drift vector components. Only for these an equivalent cubic polynomial expansion will be specified. For the constitutive equation of the equivalent system the following cubic expansion is applied, fulfilling the same symmetry conditions as the original system

d~.eq(X ) = g()c, z) = k(3c, z ) f c = blfc + b2z q-- dl.~ 3 + d2x2z .-{-d3xz 2 + d4 Z3

(20)

375 Note that due to the symmetry condition (6) the constant term and the quadratic terms in (18) are identical zero. The expansion coefficients in (20) form the tensor components B32, B33, D3222, D3223, D3233, D3333. Forming the error c = d ~ - d~,eq the least mean square conditions (8/8b 0 E[~ 2] =0, (8/8b2) E[c 2] = 0, etc. result in the following system of linear equations for the determination of the expansion coefficients ~20

~11 /~o2

E[£g]

~40

/'t31

/~22

~13

~31

/~22

~13

~04

bl b2

~60

~51

~42

~33

dl

E[zgl E[£3g]

]'L42

~33

~24

d2

EI;zg]

/~24

l'£15

d3

E[gzEg]

/x06

d4

E[ z3g]

symm.

(21)

where ~ j = E[(.~)'(z) j] and g(.~, z ) = k()~, z)&. For the fourth equation the following expansion is used d~,eq=g(£, z ) = h ( £ ,

z)£-E[h(:~,

z)£]

= a 1 + b1.¢c+ b2z + c1:~2 + c2.~z + c3 z2 + dl£ 3 + d2yc2z + d3.~z 2 .-F d4 z3

(22)

From the asymmetry condition (7) it follows that b a = b 2 = d I = d 2 = d 3 -- d4 - 0. The remaining four coefficients, forming the tensor components A4, C422, C423, C433 are then determined from the following system of linear equations 1

/X2o

/~

/~o2 l [ a l ~31 ~t22 ] [ Cl /LEE /~13 C2 =

~t40

symm.

/t04

ca

0

E['x2g] E[Yczg]

(23)

E[z2g]

The remaining tensor components of Bim are equal to the linear components of the matrix A as given by (2b). The remaining components of A i, Cimn and Dimnp in (19) are all equal to zero. /xij, i + j ~ 4 in (21) and (23) are obtained from the differential equations (19). The expectations on the right hand sides, as well as /~ij, i + j = 6 are obtained using an approximate pdf f ~ ( £ , z). This is assumed in the form of a two-dimensional Gram-Charlier series with a Minai-Suzuki modification [1,2] f~(£,z)=f°(2,

z)+3(z-1)f

+3(-z 1

fxz( , z ) = - oooo

f°(Yc, u) du

o - 1)f1-1 fx~(.;c,u) du,

E c.I-I, i+j=0

(~, z ) ~ R

Hj

× [ - 1 , 1]

(24)

(25)

376 where [i/21

Hi(x)=

tx

i

~., ( - 1 ) ai,,,x -2",

at'"

i !2 - '~ a!(i-2t~)!

(26)

or=0

1

z

(27)

H~(.) is the Hermite polynomial of the ith order. The expectation E[-]0 in (27) is with respect to the pdf f ° ( 2 , z). As a consequence of (6) the joint pdf f ° ( 2 , z) must fulfill the symmetry condition f°z(2, z) = f o ( _ 2, - z). From (27) it then follows that c°j = 0, i + j odd. Further, c ° = 1, c0°2= coo = 0. The free parameters of the tentative pdf (24) are then o°, o° and the non-trivial expansion coefficients c°1, c° , c31, o c°2, c13, o cO. These are determined so that (25) provides the same covariances and joint central moments of the 4th order for 2 and z. In order to prevent negative side loops of the approximate joint pdf c°4 is assumed equal to zero. Analytical results for/.t~j, i + j = 6 can be derived whereas the expectations on the right hand sides of (21) and (23) can be reduced to one-dimensional quadratures [3]. Equivalent linearization corresponds to Ci,,, . = Oimnp 0 in (19). This means that b 1 and b 2 are the only non-zero coefficients in (20), and all the coefficients in (22) are zero. Consequently, damage indicators for which the symmetry condition (7) is fulfilled cannot be analysed by equivalent linearization techniques. =

5. CUMULANT-NEGLECT CLOSURE TECHNIQUES The system of moment equations (19) is not closed because joint central moments ~mnijk and I~mnpijk of the 5th and 6th order, which enter the right hand sides of (19b) and (19c), are not provided by the system. Consequently, ~ijkZ,, and I.£ijklmn must be expressed by the joint central moments available. One approach is to assume that joint cumulants above the 4th order are negligible. This results in the following relationships for the joint central moments of the 5th and 6th order [10] = lO{

}s

~ijk,,.. = 15(/Zifl~k,m. }s + 10{/Zijk~,,.. }s -- 30{ ~iflZk,/~,.. }~

(28) (29)

The symbol {-}s indicates a symmetry operation producing the arithmetic mean of all terms similar to those indicates, obtained by permuting all free indices. In the following, (28) and (29) will be referred to as the conventional cumulant-neglect closure scheme. This approach is appropriate for joint probability densities which do not deviate too much from normality, or at least are of the continuous type. However, under heavy yielding the finite probability of yielding, formally appearing as delta functions in the pdf of the z-component, implies that significant joint cumulants of arbitrary order must be present. A modification of the cumulant-neglect closure scheme considering this problem is suggested.

377

Introduce the state vector X 0 = (x, Xp, u, t~). The following joint pdf is considered

fx(X)

=f°(x)+

f°o(Xo)B(z-

1)~°°f°(£, u)du -1

0

+/°o(Xo)8(-z- 1) f_ff;:(x, u) du,

[-1, 11

(30)

f°(x)

is a continuous auxiliary pdf and f°o(Xo) is the corresponding marginal joint pdf of X 0. From (30) the following expression can be derived for the joint central moments of the order

L=i+j+k+l+m+n

E[xiycjT"k(xp) `umftn] = E[x'YcJzk(x;)tu"ft"]o+ E[xi(x;)lumftn]o" ejk

(31)

where ej,=

55

2 J ( 1 - zk)f~:(£,° z) d 2 dz +

S_=S_-?

( ( - 1 ) I` -

z~)f°(Y¢,

z) d £ dz

(32)

and Xp = Xp -/lxp. The joint central moments E[xiYcJzk(Xp)tum~x"]oare, with respect to the pdf, Since this pdf is of a continuous type these expectations for L = 5, 6 can all be expressed by joint central moments with respect to fO(x) of the 2nd, 3rd and 4th order by the conventional cumulant-neglect closure scheme (28), (29). Finally, joint central moments of the 2nd, 3rd and 4th order with respect to f°(x) can be related to corresponding joint central moments as calculated by (19) by means of (31) for L = 2, 3, 4. This approach will be referred to as the modified cumulant-neglect closure scheme in the following.

f°(x).

6. NUMERICAL EXAMPLE The following data are assumed a~ = 1.0,

~ = 0.01,

f~= 2.0, b = 0.2,

7=0.50, TO= 9.0,

a

= 0.0,

a 1 = 0.0,

a 2 = 1.0 (33)

y = 2.0

Notice that ~, specifies the standard deviation of the displacement of the corresponding linear oscillator in proportion to the elastic limit displacement. Consequently ~, = 2 indicates a strong excitation, a = 0 represents the special case of an ideal elasto-plastic oscillator. To check the validity of the approximations response moments have been obtained from Monte Carlo simulation based on averaging over 5000 independent sample curves, each obtained by numerical integration of (2). The step length of the time integration has been selected as T/50, where T = 2~r/o~ is the linear u n d a m p e d eigen vibration period of the oscillator. Investigations have shown that performing more than 1000 sample curves tends to stabilize the response estimates. In Figs. 2-6 the time-dependent variation of the mean of Xp and the standard deviations of x, ~, z, Xp have been shown. Curve a corresponds to simulation, curve b to conventional cumulant-neglect closure, and curve c corresponds to the modified cumulant-neglect closure scheme. Time has been normalized with respect to T.

378 o~

0 x

1.0 1.2

a

c

0.8

¢/

0.6

0.8

°o

o

0.4 0.4 0.2

t/r 0.0 0

I

I

I

I

I

~

2

4

6

8

10

t/T 0.0 2

Fig. 2. S t a n d a r d deviation ox.

4

6

8

l0

Fig. 3. S t a n d a r d deviation ox.

The general conclusion that can be drawn from these results is that the modified cumulant-neglect closure scheme provides more accurate predictions than the conventional cumulant-neglect closure scheme, assuming the simulation results to be sufficiently accurate. From the physics of

1.0 0.8

t

°xp

Oz

2.0

1.5

~

0.6 1.0 0.4

~,,../..-

0.5

0.2

"/"'

t/T 0.0

I

I

I

2

4

6

I

i

8

1

o°°f.

0.0

10

t/T

I

I

I

I

I

2

4

6

8

10

Fig. 5. S t a n d a r d deviation axp.

Fig. 4. S t a n d a r d deviation %.

IdXp 3.0 2.5 2.0 1.5 1.0 0.5

L

0.0 2 4 6 Fig. 6. M e a n value g ~ .

8

10

379

the problem and the inherent properties of the modification (30) this is believed to be the case for any system undergoing heavy yielding.

7. CONCLUSIONS An analytical method for hysteretic structures based on an equivalent cubic expansion of the constitutive equation and the damage indicator equation has been formulated. The system of joint moment equations is closed by a conventional cumulant-neglect closure scheme and by a modified version, taking into consideration the finite probability of yielding states. Based on an analysis of a single-degree-of-freedom bilinear oscillator subjected to filtered white noise, and applying the non-dimensional accumulated numerical plastic deformation as damage indicator, it is concluded that both closure schemes provide acceptable results for the time dependent lower order moments of the damage indicator. However, the modified cumulant-neglect closure scheme is superior to the conventional version when heavy yielding takes place.

8. REFERENCES 1 R. Minai and Y. Suzuki, Seismic reliability analysis of building structures, Proc. ROC-Japan Joint Seminar on Multiple Hazards Mitigation, National Taiwan University, ROC, 1985, pp. 193-208. 2 Y. Suzuki and R. Minai, Seismic reliability analysis of hysteretic structures based on stochastic differential equations, Proc. 4th Int. Conf. on Structural Safety and Reliability, Vol. 2, 1985, pp. 177-186. 3 K.J. Mork, Stochastic response and reliability analysis of hysteretic structures, Structural Reliability Theory, Paper No. 56, University of Aalborg, Denmark, 1989 (Ph.D. thesis). 4 Y.K. Wen, Method for random vibration of hysteretic systems, J. Eng. Mech., ASCE, 102 (1976) 249-263. 5 T.K. Caughey, On the response of non-linear oscillators to stochastic excitation, Probab. Eng. Mech., 1(1) (1986) 2-4. 6 P.-T.D. Spanos, Formulation of stochastic linearization for symmetric or asymmetric M.D.O.F. nonlinear systems, J. Appl. Mech., ASME, 47(1) (1980) 209-211. 7 T.T. Baber, Nonzero mean random vibration of hysteretic systems, J. Eng. Mech., ASCE, 110(7) (1984) 1036-1049. 8 S.R.K. Nielsen, K.J. Mork and P. Thoft-Christensen, Response of hysteretic multi-storey frames under earthquake excitation, Earthquake Eng. Struct. Dyn., 18 (1989) 655-666. 9 L. Arnold, Stochastic Differential Equations: Theory and Applications, Wiley, New York, 1974. 10 R.L. Stratonovich, Topics in the Theory of Random Noise, Volume 1, Gordon and Breach, London, 1963.