Rational study of hysteretic systems under stationary random excitation

RATIONAL STUDY OF HYSTERETIC SYSTEMS UNDER STATIONARY RANDOM EXCITATION FAROUK BADRAKHAN Collegeof Engineering and Petroleum. Kuwait University, Kuwait

“trace” method together with a new definition of the normalization factor of the response d~st~bution density related to a probable maximum peak value of the response, lead to quasi-deterministic resufts which seem to be more reatistic and more rational than thosegiven by other known methods. in the particular case of bilinear hysteresis, it is shown that the statistical properties obtained are bounded by those of two limit linear cases. Abstrmt-The

1. INTRODUCTION

Structures under dynamic loading usually exhibit hysteretic behavior especially when the material is stressed beyond its elastic limit or when joint slippage or relative sliding between different components occur. In these cases hysteretic effects are the dominant factors in limiting the vibration amplitudes of structures. The response of a hysteretic system to stationary random excitation has been extens~veiy studied during recent years [I - 153. The purpose of the present paper is to apply the “trace” method introduced in f t2] in order to obtain characteristic results concerning the response of a system with hysteresis to a stationary white random excitation having a Gaussian distribution. The “trace” method is applicable to any hysteretic shape, but it is used here for the special case of a bilinear hysteresis. The method preserves the non-linearity of the system stiffness. Only the hysteretic damping is approximated and an equivalent viscous damping, function of amplitude, is defined assuming that the response is a narrow-band process (sinusoidal with a modulated amplitude). The direct writing of a Fokker-Planck equation for the density function of displacement and velocity is justified by satisfying the conditions stated by Caughey [17], and distribution densities of the response are obtained in a simple analytical closed form. The study conducted in this paper is different from those done previously in the fact that general results concerning distribution densities are exploited to define a probable maximum peak value of the motion. The statistical properties corresponding to this probable maximum are characteristic of the hysteresis shape and of the level of excitation. In this manner, the study of random vibrations will be a quasideterministic one. When applied to the case of bilinear hysteresis, the results show that the statistical characteristics of the response are always bounded by those of two limit linear cases. A Monte Carlo simulation permits to check the accuracy of the method. 2. GENERAL

For a single degree-of-freedom equation of motion is of the form

CONSIDERATIONS

system with hysteresis, under random

J f m.J3 f _f(Y> $1= sC&

forcing, the

m

where C is the viscous damping ratio, o is the natural frequency, &) is the random excitation, supposed to be a white random process with a Gaussian distribution, andf(y, j) is the hysteresis function. f&j) is a function of the displacement y and of the sign of velocity Jo(Fig. 1) and it may be written [I23 in the form

f(u, 9) = dy) + NJ) sgn3,

(2)

or in the form

S=f*. 3t5

(3)

F. BADRAKHAN

316

Fig. I. Bilinear hysteresis (I; k).

with f+ =

AY)+

if j > 0,

NY)

if Y < 0.

f- = gtY)S h(Y)

If the hysteresis cycle is centered, the equation F = g(y) = 1/2(j+ + f-) is the equation of the “trace” of cycle. g(y) is an odd function of y, while h(y) = 1/2(f+ -f_) is an even function. According to the “trace” method [12], the cycle defined by F = h(y)sgn j is approximated by an ellipse, of the same amplitude and enclosing the same area, representing an equivalent viscous damping. Due to this approximation, equation (1) can be rewritten in the simpler form j; +

w + g(Y)

=

q(r).

(4)

with B = 215 + B,,

(5)

where B, is the factor of the equivalent viscous damping. It can be easily evaluated from the cycle area through the relation

s A

I~B,wA’

= area of cycle = 4

0

h(Y)dY,

(6)

where A is the response amplitude. For the particular case of a bilinear hysteresis (l;k), illustrated in Fig. 1, B, is given by the expression

B, =

4(1 - k)(A - 1) mA*

*

(7)

Expressions (6) and (7) are evidently valid for a steady sinusoidal motion of the system. Under a white random excitation (broad-band process), with the assumption of a moderate non-linearity, the response of the hysteretic system is supposed to be a narrow-band process (sinusoidal with a modulated amplitude). It is assumed also that the amplitude modulation is sufficiently slow that equations (6) and (7) remain valid for each cycle of the random response of the system. The frequency o and a characteristic amplitude to be used in these equations are defined below. For a bilinear hysteresis (l;k), the function g(y) represented in the force-displacement

Hysteretic

317

systems

plane by the “trace” is defined as follows:

for 0 < y < 2 - A, g = l/2(1 + k)y - l/2( 1 - k)(A - 2)

A<2

g = I

for 2 - A < y < A;

63)

for 0 G y G A - 2,

ky

g = f/2(1 + k)y - l/2(1 - k)(A - 2)

for A - 2 G Y G A.

In all other hysteresis models considered so far [6, 13, 14, 163, the function g(y) can be easily defined; it is an odd function derivable from a potential. Furthermore, g(y) may be replaced by a polynomial obtained by appropriate fitting or from the first terms of Maclaurin’s expansion in the form: g(y) = a,y + a,y3 + *** . It is to be noted here that the concept of writing the equation of motion for a hysteretic system in the form (4) seems to be similar to that of [S]. Actually, it is completely different since the non-hysteretic function in the present case is not an arbitrarily selected function, but it is the well defined trace of cycle. Since g(y) is an odd function and the energy dissipation is due to a damping force proportional to velocity, a Fokker-Planck equation may be written [17]. If only the stationary response is sought, this equation is of the form

ap -$a (CBu + &?tY)lPI+

“ay=

.@

where p = p(y,u) is the joint distribution density of displacement and velocity u = 3, S is the power spectral density (PSD) of the excitation q(t) defined from E[q(t)q(t - r)] = 2nS6(r), where EC...] is an ensemble average and 6(r) is Dirac’s impulse. Since q(t) is considered here as a white noise, S is constant for a given level of excitation. The unique solution of equation (9) is of the form (see [17])

P(Y,

4 = P(Y)*~(4 = C exp

(10)

where

G(Y) =

&)dx

(11)

is the potential energy, and 7LS 0,’ = B

(12)

is the velocity variance. The definition of C: by equation (12) is identical to that obtained by Karnopp’s power balance method [ 181. The mean square value of displacement is obtained from the relation

CT2 02’

(9 = 2

Y

(13)

c,’ is related to the amplitude A of the response, and for bilinear hysteresis (1; k), for

318

F.

BADRAKHAN

instance, it is defined by 4(i --&(A - If nwA2

The constant C in equation (10) is a normalization

1 ’

factor. It is defined by

In this definition of C, it is understood that u and y can take, during motion, all possible values irrespective of the excitation level. Actually, this is not the case, especially for the displacement y. In order to study the inRuence of the excitation level (defined by the power spectral density S), the normalization constant C wilf be replaced by another approximate one C* given by the relation

where A* is a probable maximum peak value of y defined by Frob~bility (peak > A*) c E.

(171

The number E can be chosen as small as desirable. The value E = 0.001 will be adopted in the subsequent applications in this paper. The importance and originality of the new definition given to the normalization factor C* from equations (15) and (17) will be demonstrated by the results obtained which are very characteristic of the hysteresis systems under investigation. Since 6, is a function of A, and in order to make a useful comparative study, it is obvious that one should consider the statistical properties of the response for a particular well defined value of A. Generally, in a major class of first-passage problems, a critical peak value A, is defined, and subsequent analysis of the response is performed in terms of this critical amplitude A,; but if the problem of random response in its general sense is under investigation, it is logical to consider the probable maximum peak A* defined from equation (17) as characteristic and intrinsic for the system under a given excitation level. Vice versa, A* is an identi~er of the excitation level for a given system. For this reason, A* is used in what follows in order to study the influence of both hysteretic and viscous damping and of the effect of the PSD of excitation on the response. A particular cycle of amptitude A* will be considered as a characteristic cycle corresponding to the given level of excitation S. This will make the study of random vibrations quasi-deterministic. Setting a,(A*) =

(18)

6,

and since p{y) and p(v) are separabIe, one can write

p(u) = 1exp

I

aJ2n

p(y)=$exp

i

-F

-$ 1 , I

,

(1%

WI

319

Hysteretic systems

with

(21) The probability function of equation (17), used to define A*, will be established in the following section. The factor Y defined by equation (21) can be considered as a function of the PSD of excitation: Y = Y(S) and, consequently, an intrinsic distribution of displacement is defined for every level S. For a bilinear hysteresis (1; k) the expressions of 9 are given in Appendix A; those for a polynomial trace (g = a,y + a,y3 + **a)are given in Appendix B. 3. STATfSTICAL

PROPERTIES

OF THE

RESPONSE

The statistical properties of the response are described in [!23. They are reconsidered here in forms adaptable to applications. The response process crosses, with a positive slope, the level y = A a number of times per second. The expected value N: of this number of times is given by Rice’s formula [19-J:

s 33

N:

(22)

MY~I,=,U~G

=

0

which yields, after integration,

(23) For A = 0, the usual pseudo-frequency

is obtained (24)

Since G(0) = 0, equation (2) may be rewritten in the form

G

=Niexp

(25)

The circular frequency w to be used in equations (5)-(7) is taken as o=wg+

=2xN;.

(26)

The probability that the amplitude of the response process is larger than a given value A may be written in the form (reference [20))

P = P(peak > A) = f$

= exp

,

(27)

0

taking account of equation (25). The peak distribution density is given by

~&f)=~= G(A) 1 dP

d -Gexp

[

-a2

.

c3-0

Equation (27) is used to determine the probable maximum peak value A* mentioned above

320

F. BADRAKHAN

in equations (16)-(21). Specifically, A* is defined from the equation

P(peak > A*) = exp

= 0.001.

Since d is also a function of A*, few iterations are required to determine A* with the desired accuracy. The envelope distribution density of response may be determined using Crandall’s approach [20] which gives, for this density,

P,(A)= - N,+T(A)

-&{““P[ - y]},

(30)

where T(A) is taken as the period of the free undamped system J + g(y) = 0. Defined in this manner, the envelope distribution density is similar to the energy-based envelope density and is normalized, [21,22]. Expressions of the function G(A) for bilinear hysteresis (1; k) are given in Appendix C. 4. APPLICATION

TO THE

CASE

OF

BILINEAR

HYSTERESIS

AND

DISCUSSION

Figures 2-8 show applications of the above general results to the case of bilinear hysteresis (1; k). In Fig. 2, variations of the frequency Ni are represented, in terms of the yielding ratio A for various values of k. Curves representing A* in terms of k, for various values of the viscous damping ratio and for a particular value of the level S of excitation are shown in Fig. 3. Figures 4 and 5 represent variations of the reduced distribution density of y, ap(y), in terms of the reduced displacement :. It is shown that these curves tend to the normal (Gaussian) distribution defined by

OP(Y) =

&exp[--li21 ,

(31)

when k tends to one, for a given value of S. In fact the curves ap(y) are always bounded by the normal distribution (31) and by another normal distribution, depending on k, defined by the density function

.

0

1

2

SlMULATlON

3 (A;

Fig. 2. Frequency Ni vs level A. Results of digital simulation are shown.

321

(k) Fig. 3. Maximum peak value defined from P(peak > A*) = 0.001. vs k.

k *

0.2 - 0.8

---- SIMULATION

(Y/W)

Fig. 4. Comparison of distribution density of displacement with the results of digital simulation.

(32) which is attained when k decreases for a given value of S or when S increases for a given k. The two limit Gaussian curves are represented in Fig. 5 for the case k = 0.25, where it is shown that when the excitation level changes, the distribution curves of Hy) are bounded by the curve defined from equation (31) and that defined from spot) =

GFk51)

. The latter is reached for a reiatively low level of excitation (S > 0.05).

Th& result concerning the existence of two limit distributions is predictable because of the bilinear character of the system and since it is known that the distribution density of the response of a linear system is Gaussian. But, this result is obtained here for the first time because of the new definition given to the normalization factor and because of the approach used which preserves the non-linearity of the system while all previous studies [l, $7, lo] are based either on a linearization approach or on averaging techniques which mask the bilinear character of the system. The results of a Monte Carlo simulation are shown in

322

F.

BADRAKHAN

ke

0.25

Y t 0.00

s:

0.01 - 0.10

:.:.‘:,E SIM‘,,_ATIO,., P

b

a2

a0 1

0

2

L (Y&7

5

?

Fig. 5. Distribution density of displacement for k E 0.25 and two levels of excitation S = 0.01 and

1. 1

S = 0.10. Convergence

ady)

toward -Y2

[ 8oz

= Gexp

the two Gaussian

distributions

q(y) = --!--exp -5 &c [

1

or

and the results of the Monte Carlo simulation are shown.

Figs 3, 4 and 5. For the simulation of random white noise, normally distributed random numbers are generated at equal intervals of the dimensionless time r = or. The interval As is chosen equal to $

(as-in [7]). Random numbers Ni are assigned at these intervals

cNi

c is a normalization factor depending on the power spectral in the form q(ri) = -where 6 density S of excitation. This approach adopted for the digital simulation of white noise is similar to that of [23]. The ensemble size taken for each case is 2000. With this discrete simulation of the white noise, the original equation of motion (if is then integrated by steps using the Runge-Kutta method, The figures show that the approximate results using the trace method are in good agreement with those obtained by digital simulation. Figures 6 and 7 represent variations of a&4), reduced distribution density of peaks, in terms of the reduced amplitude A/a for various values of k, S and {. All these curves are bounded

0.6

k x 0.2-0.3-O.& 0.5 -0.6-O-7

ai

A2 1

Fig. 6. Distribution density of peaks without viscous damping (; = 0). Convergence toward the A is shown. Rayleigh distribution up, = ; exp - 9 i

Hysteretic

323

systems

0.6

_

0.‘

a? e

0.3

0.2

0.0

0

I

5

2

6

Fig. 7. Distribution density of peaks for k = 0.25 and different toward the two Rayleigh distributions oP = A/aexp[ -AL/Z&j shown.

7

8

levels of excitation. or up,, = A/4oexp

Convergence [ - A*/Va*] is

by those of the two limit Rayleigh distributions:

=aexp --2oz 1, A2

WP

[

(33)

and ap, = kaexp

[

-g

. 1

(34)

The first one (33) is reached when k increases toward one, for a given S, or when S decreases for a given k, while the second limit Rayleigh dist~bution (34) is attained when k decreases for a given S or when S increases for a given k. It is shown in Fig. 6 how the first limit distribution is approached when k increases. The second limit dist~bution, function of k, is not shown in this figure for the sake of clarity. In Fig. 7, the two limit Rayieigh distributions are represented for the particular value k = 0.25; the second limit is reached for S > 0.10. It is again worthy to note that the existence of these two limit Rayleigh distributions is a direct consequence of the expression adopted for the normalization factor and of the bilinearity of the system. This is predictable since it is known that the distribution density of peaks for the response of a linear system is of the Rayleigh type. This result is obtained here for the first time, however, because of the direct approach used which preserves the non-linearity of the trace of the hysteresis cycle. The distribution density of the response envelope is very close to that of peak values as shown in Fig. 8, where results of digital simulation for the peak values are also shown. This means that the distribution curves of the envelope are also bounded by the two Iimit Rayleigh distributions defined by equations (33) and (34). 5. CONCLUSION

The method used here to study the response process of a hysteretic system is based on the concept of the “trace”, defining the shape of the hysteresis cycle, and on the Karnopp’s definition of the mean square velocity from the power balance concept. It takes into account the non-linearity of the system and for a bilinear hysteretic system, the results are bounded by those of two limit linear systems. For the distribution density of displacement, the two limits are Gaussian distributions, while for the distribution density of peaks or envelope, the two Iimit results are Rayleigh distributions. The results are obtained in closed forms, are relatively simple and can be applied to any shape of hysteresis. In the worst case, integrations may have to be performed numerically. NC,” z?i“D

324

Fig. 8. comparison of distribution density of peaks and of the response envelope for S = 0.50, < = 0.01 and two values of the relative yielding slope. Results of digital simulation are also shown.

Acknowledge~eats-This research was supported by the University Research Council (Grant EM 19). Computer facilities were obtained from the Kuwait University Computer,Services. The author wishes to thank engineer Ghassan Shaban for his assistance in implemenfing and running the computer programs,

f. T. K. Caughey, Random excitation of a system with biiinear hysteresis. J. appt. Meek 27,649~652 (1960). 2. W. D. fwan and L. D.-Lutes, The response of a bilinear hysteretic system to stationary random excitation. J. Acous. Sot. Am. 43, 545-552 (1968). 3. D. Karnopp and R. Brown, Random vibration of multidegr~-of-fr~dom hysteretic structures. J. Acous, Sot. Am. 42,54-59 (1967). 4. H. Takemiya and L. D. Lutes, Stationary random vibration of hysteretic systems. J, Eng, Me& I)%., ASCE 103,673-687 (1977). 5. L. D. Lutes, Approximate technique for treating random vibration of hysteretic systems. J. Acoust. See. Am. 4% 299-306 (1970). 6. Y. K. Wen, Method for random vibration of hysteretic systems, J. Eng. Meek. L)io., AXE. 102, 249-263 (i976). 7. 3. B. Roberts, The response of an oscillator with bjlinear hysteresis to stationa~ random excitation, J. uppl. Me& 45923-928 ( 1978). 8. R. L. Stratonovitch, Topics in the Theory of ~u~dom Noise, Vol. 2. Gordon & Breach, New York (1964). 9. R. 2. Kbasminskii, A limit theorem for solutions of differential equations with a random rjght-band side, Theor. Prob. Applic. I f , 390-406 (I 966). 10. W. D. iwan and P. T. Spanos, Response envelope statistics for non-linear oscillators with random excitation, .i, app~. Mech. 45, 170-174 (1978). t 1. P-T. D. Spanos, Hysteretic structural vibration under random load. J. Acous. SOC.An 65,404-410 (1979). 12. F. Badrakhan, Etude approch&e d’un sysdme B hyst&%s. Int. 3. eon-li~eur Me&. 12, i-12 (1977). 13. R. Bout, Forced vibration of mechan~~i systems with hysteresis, Pr~eedjngs of 4th ~on~~e~~e on NonLineur Oscillation, Prague, Czechoslovakia (1967). 14. N. Edelmann, Model for calculating magnetic hysteresis 1oops. Siemens Forsch. ~ac~~ckt~~gs~er~ckte 5, 123128 (1976). 15. V. V. Bolotin, Tire Dynamic Sru~i~jr~ of Ei~ic S~@e?ns.~otden-Day, San Francisco (1964). 16. F. Badrakhan, NouvelIe methode d’etudc des systemes d hyst&rcisis,Th&se de bcteur b-Sciences, 3esanCon University” France (1972). 17. T. K, Caugbey, Derivation and appljcatiun of the ~okker-PIanck ~uation to discrete non-liner d~~arn~c systems subjected to white random excitation. J. Acaru. Sac. Am. 35, 1683-i692 (1963). 18. D. Karnopp, Power balance method for non-liner random vibradon, J. appt. Me& 34,212-214 ( .,. . _ _I ., _L __?__,.__1..2__d+ ___,a__--:^A i- P-t.%..._,4 II,...“_” ,... hf.rirn “..,? ~mrI.n‘.,:,.D

_. -. _

,

pp. t33-194. Dover, New Yor< (1954). 2% S. H. Crandail, Zero crossings, peaks and other statistical measures of random responses. J. Acous. See, Am. tf.

35, 1693-1699 (3963). S. I-L Cfandall, Random forcing of non-linear systems. In Les Vibrntions For&es dans le.5 Systems NonLineaires. pp. 57-68. CNRS, Paris (1965).

22. T. K. Caugbey, On the response of a class of non-linear oscillators to stochastic ex~tation. In &es ~b~ur~o~s For&es dans les Systems Non-Lineuires, pp. 393-405. CNRS, Paris (1965).

23. P-T. D. Spanos, Monte Carlo simulations of responses of non-symmet~c exultations. Comput. Struct. 13, 371-376 (198 I).

dynamic system to random

325

Hysteretic systems APPENDIX

A

Factor of normalization Y for bilinear hysteresis (1;k): A d I -. G(y) = 4~‘; Y = a$ierf>.

-f f
G(Y) = iY2 G(y) = f(l + k)y* - f(i - kXA - 2)~ + f(k - 1x2 - A)’

4p = u$iierf(%)

for 2 - A d y c A;

+ 2o&{erf[‘~J~]

- erf[~]~exp~’

A,2+

for 0 Q y C 2 - A,

t:l’,:“‘]. for 0 6 y c .4 - 2,

G(Y) = tky* G(y)=~(l+k)y2-~(l-k)(A-2)y+~(l-k)(A-2)z

forA-2CyGA;

In these expressions A* is replaced by A for convenience and erfx = %_oe f

-r’dt = error function.

APPENDIX

B

Factor of normalization Y for a polynomial trace: i ati_ ty2i-1. 1-1

g(y)=a,y+a,r3+...+a~=

with

CL’=

c

2b2,_, A”- ‘,

i=i b._

28

x021-I 6’ ’

*

The series expressing .Y is convergent in virtue of d’Alemb&t’s rule:

APPENDIX

C

Expressions of G(A) for bilinear hysteresis (1; k): A<

l-+G(A)=fA’;

1 Q A *I 2 + G(A) = f(2k - 1)A’ + 2(1 - k)A f k - f; A 2 2 -. G(A) = +kA2 f

1 - k.