Statistical random response analysis and reliability design of structure system with non-linearity

ARTICLE IN PRESS

Mechanical Systems and Signal Processing Mechanical Systems and Signal Processing 19 (2005) 1135–1151 www.elsevier.com/locate/jnlabr/ymssp

Statistical random response analysis and reliability design of structure system with non-linearity Byungyoung Moon, Choon-Tae Leea, Beom-Soo Kangb, Byeong Soo Kimc a

Department of Mechanical and Intelligent Systems Engineering, Busan National University, 30 Changjeon-dong, Keumjeong-ku, Busan 609-735, Korea b Department of Aerospace Engineering, Pusan National University, Kumjung-ku, Pusan 609-735, Republic of Korea c School of Mechanical and Automotive Engineering, Inge University, 607 Obang-dong, Kimhae, Kyungnam 621-749, Korea Received 21 May 2001; received in revised form 6 October 2003; accepted 25 May 2004 Available online 23 November 2004

Abstract In this paper, a method to analyse a complex non-linear structure system under random excitation is proposed. First, the actual random excitation, such as earthquake, is approximated to the corresponding Gaussian process for the statistical analysis. When modelling the overall non-linear system, the overall system is divided into substructures and these are approximately transformed into modal coordinates with the random excitation. The modal equations of the overall system are expanded sequentially according to the non-linear order. Then, the perturbed equations are synthesised into the overall system and solved in a probabilistic way. Several statistical properties of a random process that are of interest in random vibration applications are reviewed in accordance with the non-linear stochastic problem. The obtained statistical properties of the non-linear random vibration are evaluated in each substructure. Comparing with the results of the numerical simulation proved the efficiency of the proposed method. r 2004 Elsevier Ltd. All rights reserved.

1. Introduction In recent years, the trend in mechanical systems has been toward high-speed and lightweight ones in many industrial machines. These conditions can cause trouble of a non-linear vibration in Corresponding author.

E-mail address: [email protected] (B. Moon). 0888-3270/$ - see front matter r 2004 Elsevier Ltd. All rights reserved. doi:10.1016/j.ymssp.2004.05.003

ARTICLE IN PRESS 1136

B. Moon et al. / Mechanical Systems and Signal Processing 19 (2005) 1135–1151

mechanical systems. Hence, it has become important to consider the non-linear characteristics in vibration analysis, such as design of the structure system. For a non-linear system, exact solutions are generally not possible and approximate solutions can be obtained numerically. This method can be extremely time consuming and is not practical, especially for large degree-of-freedom (DOF) system. Iwatsubo et al. [1] have proposed a new method to analyse the vibration of a non-linear mechanical system. Moon et al. [2,3] have reported a study on the vibration of the mechanical system to analyse the dynamic problems of a non-linear multi-degrees-of-freedom (MDOF) systems. They developed the substructure synthesis method (SSM) technique to reduce the overall size of the problem for the non-linear structure, and obtained approximate solutions of the nonlinear system using an approximation method. On the other hand, it is necessary that a high-speed system used for the jet engine of an aircraft, power plant turbine and ship propulsion system, promptly pass a critical speed. Accordingly, the casing is often modelled elastically to decrease the critical speed. When random excitation excites such a mechanical system, the casing is excited to contact with the rotor and there is a danger that the bearing will be damaged. Therefore, the investigation of the random response of rotating machinery is very important from the viewpoint of disaster protection. Soni et al. [4] and Srinivasan et al. [5] have reported the random analysis of a rotor system using the response spectrum method and time response method in the deterministic system. Rao [6] has reported the method to analyse the earthquake response of the high-speed rotating machinery. They used real earthquake data to analyse the linear response. From the viewpoint of the dynamic response of the mechanical system against random excitation, it can be treated as a stochastic problem. There are analytical studies about the stochastic system dealing with single-degree-of-freedom (SDOF). Analytical solution is possible to find the power spectrum density (PSD) of the linear system, subjected to narrowband excitation by the method of measuring filters. [7] Dimentberg et al. [8] obtained an exact solution for the spectral density of the response of a certain specific non-linear SDOF system, excited by a white noise force. Some researchers have studied the analytical method to estimate the PSD response with nonlinear restoring force. Bouc et al. [9] presented an analytical procedure to estimate the PSD response of a weakly damped oscillator with a non-linear asymmetrical force under external stochastic wide-band excitation. Bouc et al. [10] proposed an analytical method for approximating second-order statistics (mean, covariance matrix, PSD matrix function) of response of the MDOF mechanical systems with non-linear strongly restoring forces and small damping, when the excitation is known only by its PSD matrix. Bouc et al. [11] presented a method for estimating the PSD matrix of the stationary response of lightly damped randomly excited MDOF mechanical systems with strong non-linear asymmetrical restoring forces. However, the vibration of non-linear rotor systems, which are modelled by MDOF, under a random excitation force has been less investigated so far. Moreover, the vibration analysis of a non-linear rotor–bearing–casing system utilizing a statistical approach to a seismic wave is not found in past research. Therefore, this paper proposes a numerical method for non-linear vibration of a mechanical system against a random excitation by applying the statistical method. The actual random excitation, while regarding earthquake excitation as a stationary random process, is approximated

ARTICLE IN PRESS B. Moon et al. / Mechanical Systems and Signal Processing 19 (2005) 1135–1151

1137

to the corresponding Gaussian process. Then, several statistical properties that are of interest in non-linear random vibration applications are reviewed in accordance with the stochastic problem.

2. Method of analysis The mechanical structure system with the non-linear restoring force (or moment) of the system, which is excited by earthquake, is considered as shown in Fig. 1. The excitation is regarded as a random process; hence, it is extremely difficult to obtain an exact solution. Thus, solutions can be obtained approximately. It should be noted that the solution itself for random inputs is not the ultimate goal in stochastic analysis of a non-linear system; instead, more relevant information is the statistical properties of the amplitude of the response. To elaborate, a non-linear system with the inputs, which are assumed to be a Gaussian random process, is considered. Because of the non-linear characteristic, the output is no longer a Gaussian random process; hence, the statistic characteristic of its amplitude cannot be evaluated through the probability density functions (PDF). Therefore, an adequate method to evaluate the statistical properties of the response of a nonlinear structure system should be developed. For this reason, the random excitation is approximated to the Gaussian stationary process by reasonable procedures. Then, the non-linear equation of motion reinstates with the approximated Gaussian process. After that, the perturbation theory is applied to solve the non-linear equation of motion. Finally, the statistics properties of the non-linear response are obtained. 2.1. Non-linear system excited by random process For the simple explanation of methodology, the equation of motion of the arbitrary mode of the non-linear system against random excitation can be expressed as x€ þ 2Bo0 x_ þ o20 x þ
o20 x3 ¼ X€ 0 ;

Fig. 1. Non-linear Mechanical system for analysis.

(1)

ARTICLE IN PRESS B. Moon et al. / Mechanical Systems and Signal Processing 19 (2005) 1135–1151

1138

where B and o0 are the damping ratio and natural frequency, respectively.
is a small non-linear parameter. Eq. (1) has a non-linear restoring force, which can be expressed in higher order terms of the displacement,
o20 x3 : X€ 0 is random excitation, which has spectrum density SN ðOÞ: Meaning of the stochastic analysis such as Eq. (1) is to decide the statistic information of displacement x. Generally, the statistic characteristic of the random process is decided from the PDF and PSD function of the system. Accordingly, for the probabilistic analysis of the non-linear random response, PSD and PDF of excitation forces should be obtained. It is clear that ground acceleration is inherently non-stationary, treating a typical record of earthquake-induced ground acceleration [12]. However, if the principal shock duration is limited to the period corresponding to the strong-motion portion over which the peak structural response occurs, a stationary process appears to be a good approximation. Hence, this study deals with the response of a non-linear system under ground excitation of the stationary random process by considering the strong motion duration of the earthquake. Fig. 2 shows the random excitation of the system and its PSD and PDF, which is regarded as a narrowband process. Important values of statistic properties are mean (=0.072) and peak frequency (=28.27 rad/s) and maximum value of acceleration (175.89 Gal). Fig. 2(b) shows the corresponding erratic PSD and the fitted PSD function, which show good agreement. The computed PSD is obtained from the part of earthquake data during 3–18 s, which can be regarded as a stationary process from the viewpoint of the amplitude envelope of earthquake graph. The -4

x 10 20 Spectrum (m/sec 2/Hz)

Acceleration (m/s 2)

2

1

0

-1

15

Computed PSD function

10 Fitted SN( Ω) 5 0

(a)

10

20 30 40 50 Time in seconds Taft Earthquake (1952, S69E)

60

0

(b)

20

40 60 80 Frequency in rads/sec PSD function

1.5 Gaussian distribution 1 PDF

0

0.5

0 -2

(c)

-1

0 1 Acceleration (m/s 2) PDF of iv

2

Fig. 2. Random excitation force of the system and its PSD, PDF.

100

ARTICLE IN PRESS B. Moon et al. / Mechanical Systems and Signal Processing 19 (2005) 1135–1151

1139

functional form for the spectral density function of the earthquake motion can be expressed as SN ðOÞ ¼

o2g a2 O2 ðo2g

2 2

O Þ þ

4B2g o2g O2

S 20 ;

(2)

where og ; zg and S 0 are a dominant frequency, damping ration of filter and spectrum intensity of the random process, respectively. a is a maximum value of input time history. For instance, the Taft earthquake (1952) has values of Bg ¼ 0:41; og ¼ 18:75 rad=s; a ¼ 1:75 m=s2 and S0=0.0132 m2/s4/Hz. From Fig. 2(b), it is estimated that SN ðOÞ can be applicable to solve the response statistically. PDF of earthquake input is obtained from the part of earthquake during 3–18 s, which shows a Gaussian distribution, as shown in Fig. 2(c). The part of the earthquake during 3–18 s has the statistic properties as mean (=0.0062) and the maximum value of acceleration (=175.89 Gal). Thereby, the part of earthquake excitation process can be regarded as Gaussian, stationary random process with mean zero, where the PDF can be expressed as ! €2 1 X 0 (3) PðX€ 0 Þ ¼ pffiffiffiffiffiffiffiffiffiffiffiffiffi exp  2 : 2psX€ 0 2sX€ 0

According to the procedure described above, the random excitation is expressed by the Gaussian stationary process in the analysis of the non-linear system. Then, the non-linear equation of motion can be reinstated as x€ þ 2Bo0 x_ þ o20 x þ
o20 x3 ¼ W ðtÞ:

(4)

The excitation W ðtÞ is a Gaussian stationary narrowband random process with zero mean, which has the spectral density function S N ðOÞ: xðtÞ is the stationary response of a linearly damped Duffing system subjected to a Gaussian process excitation. Here, xðtÞ is not a Gaussian process, generally. However, if the system is lightly damped and if the non-linearity is small (i.e. is,
p1), then xðtÞ is still expected to be a narrowband random process, as shown in Fig. 3. The response has the mean (=0.0044). PDF of the non-linear response shows a Gaussian distribution with mean zero. Thus, the response of Eq. (4) can be analysed statistically by applying the perturbation theory. In Eq. (4), the small parameter
can be perturbed as x ¼ x0 þ
x1 : The perturbation method replaces the solution of the non-linear random vibration problem represented by Eq. (4) with the solution of a series of linear random vibration problems with the same differential form but different inputs x€ 0 þ 2Bo0 x€ 0 þ o20 x0 ¼ W ðtÞ;

(5)

x€ 1 þ 2Bo0 x_ 1 þ o20 x1 ¼ f p ðx0 Þ;

(6)

where f p ðx0 Þ ¼ 
o20 x30 : From Eqs. (5) and (6), the approximating functions x0 ; x1 can be obtained sequentially. The zeroth-order approximation x0 is the Gaussian process. Its probabilistic characteristics can be obtained by classical methods of linear random vibration. However, the characterisation of x1 is less simple because the input is now a non-Gaussian process whose mean and covariance function are not generally available in a closed form. Moreover, the

ARTICLE IN PRESS B. Moon et al. / Mechanical Systems and Signal Processing 19 (2005) 1135–1151

1140

1.5 Linear response Nonlinear response

1 0 -1 -2

(a)

Probability density function

Displacement (cm)

2

0

5

10 15 Time in seconds Time response

20

Gaussian distribution 1

0.5

0 -2

25

(b)

-1

0 1 Displacement (cm) PDF of nonlinear response

2

Probability density function

1.5 Gaussian distribution 1

0.5

0 -2

-1

0 1 Displacement (cm) PDF of linear responsz

(c)

2

Fig. 3. Non-linear response x; linear response x0 of Eq. (4) and their PDF (when
¼ 0:3; B ¼ 0:1; o ¼ 5:23).

determination of the second moment characterisation of the approximate solution x ¼ x0 þ
x1 also requires the correlation function of x0 ; x1 ; which is not readily available. In handling this problem, the objective is the determination of the stationary mean and covariance function of the first-order approximation, as given in Eq. (6). Covariance function becomes, due to its stationary characteristic, Efxðt þ tÞxðtÞg ¼ RðtÞ ¼ Efx0 ðt þ tÞx0 ðtÞg þ
Efx0 ðt þ tÞx1 ðtÞg þ
Efx1 ðt þ tÞx0 ðtÞg; where each term can be obtained from the random vibration of the linear system as Z 1Z 1 EfW ðt þ t  y1 ÞW ðt  y2 Þghðy1 Þhðy2 Þ dy1 dy2 ; Efx0 ðt þ tÞx0 ðtÞg ¼ 0

(7)

(8)

0

Z

1

Z

1

Efx0 ðt þ tÞx1 ðtÞg ¼  0

Z

1

Z

Efx1 ðt þ tÞx0 ðtÞg ¼  0

EfW ðt þ t  y1 Þf p ðt  y2 Þghðy1 Þhðy2 Þ dy1 dy2 ;

(9)

EfW ðt  y1 Þf p ðt þ t  y2 Þghðy1 Þhðy2 Þ dy1 dy2 ;

(10)

0 1

0

where h is the impulse response function corresponding to
¼ 0: The determination of the expectations in Eqs. (8)–(10) involves lengthy but straightforward calculations of expectations of polynomials in Gaussian variables. Thus, the covariance function of the non-linear response can

ARTICLE IN PRESS B. Moon et al. / Mechanical Systems and Signal Processing 19 (2005) 1135–1151

1141

be evaluated. Then, the spectral density of the non-linear response is obtained by taking the Fourier transform of Eq. (7) Sxx ðOÞ ¼ S N ðOÞjHðOÞj2 ½1  6
o20 s2x0 RefHðOÞg;

(11)

where s2x0 is the stationary variance of the linear response x0 ðtÞ: HðOÞ is the frequency response function of the linear equation between the excitation and the displacement of response. The corresponding variance can be obtained from the covariance RðtÞ of the system by letting t ¼ 0; which is the same value of the mean square response of the non-linear vibration. Since the mean response, E½xðtÞ; is zero, the variance is equal to the second moment E½xðtÞ:   Z 1 2 2 2 fRðtÞhðtÞg dt ; (12) E½x ðtÞ ¼ Rð0Þ ¼ sx0 1  6
o0 0

The mean square value of the linear response in terms of the system response function and the spectral density of the input random process can be obtained as Z 1 s2x0 ¼ S N ðOÞjHðOÞj2 dO: (13) 0

This procedure is applied to analyse the complex multi-DOF non-linear system. 2.2. Modelling of the complex mechanical system For the analysis, the system shown in Fig. 1 is considered as a mechanical system. The rotor has non-linearity with respect to its material property. For the dynamic analysis of complex systems, the SSM can be applied. The overall system is divided into three components, i.e. the rotor is the non-linear component, the casing is the linear component and the bearing is the assembling component. The coordinate system of the rotor–bearing–casing system is shown in Fig. 1. The Or  X r Y r Zr coordinate system is fixed in rotor such that the origin coincides with the center of the shaft where the X r -axis is vertically upwards, the Y r -axis is horizontal and perpendicular to the shaft and the Zr -axis is along the shaft. The Oc  X c Y c Z c coordinate system is fixed in casing. The O0  X 0 Y 0 Z 0 is an absolute coordinate system, which is fixed in basement. The U r ð¼ X r  X c ; Y r  Y c ; Z r  Z c Þ is a relative displacement between the rotor and casing. The U c ð¼ X c  X 0 ; Y c  Y 0 ; Z c  z0 Þ is a relative displacement between the casing and basement. X€ 0 is an acceleration of the earthquake input. The acceleration of gravity is ignored for simplicity. The shaft and casing components are modelled using FEM. 2.3. Modelling of non-linear component When the rotor is modelled by the finite element method (FEM), the characteristic of the nonlinear restoring force for each element is based on the relation that the stress of the element by bending moment is represented by the sum of two terms, one that varies linearly with the strain and one that varies with the third power of the strain. The non-linear stiffness term of the rotor is determined approximately. Then, the non-linear restoring force is formulated in a complex form including non-linear coupling terms of displacement. However, as the first step of the analysis, it is

ARTICLE IN PRESS B. Moon et al. / Mechanical Systems and Signal Processing 19 (2005) 1135–1151

1142

assumed that the non-linear restoring force fRe g is represented as [2,3] fRe g ¼ ½1 K e f1 X e g þ
½K Ne f1 X 3e g; f1 X e g ¼ fxr ; yrx ; yr ; yry gT ; f1 X 3e g ¼ fx3r ; y3rx ; y3r ; y3ry gT ;

(14)

where [1Ke] and [K Ne ] are the stiffness matrix and non-linear stiffness term of the element, respectively.
is a small parameter. The superscript 1 denotes the non-linear component. ð ÞT denotes the transpose of ð Þ: The external force, which acts on the rotor, is considered as the unbalance force and the earthquake acceleration. Internal force is considered because the nonlinear component can be synthesised through the internal force with the other components. By considering the boundary conditions, the equation of motion for the non-linear component can be written as [2,3] ½1 Mf1 U€ r g þ ½1 Kf1 U r g þ
½K N f1 U 3r g ¼ f1 F u ðtÞg þ f1 F E g þ f1 F b g;

(15)

where [1M] is a mass matrix. {1Fu(t)} is an external force vector caused by the unbalance of the rotor. {1Fb} is an internal force vector.{1FE} is an external force against the earthquake acceleration f1 F E g ¼ ½1 MfIgfX€ 0 g;

(16)

where {I} is a vector which shows the direction of the earthquake. In order to apply modal analysis, modal coordinate system is introduced by using the modal matrix ½1 F of the only linear system. Then, the displacement {1Ur} in the physical coordinate system can be transformed into the modal displacement approximately as [2,3] € þ ½\1 o2 f1 xg þ
½\1 kN\ f1 x3 g ¼ f1 f ðtÞg þ f1 f g þ
f1 f g; f1 xg u E b \

(17)

where f1 f u ðtÞgð¼ ½1 FT f1 F u ðtÞgÞ and
f1 f b ðtÞgð¼ ½1 FT f1 F b gÞ are the external and internal forces in modal coordinates. f1 f E gð¼ ½1 SfX€ 0 g; ½1 S is a mode-participation factor of its mode) is the external force against earthquake in modal coordinates.
½\1 kN\  is the non-linear term in modal coordinates. According to the previous section, Eq. (17) can be expressed with the approximated Gaussian stationary random process W ðtÞ as € þ ½\1 o2 f1 xg þ
½\1 kN\ f1 x3 g ¼ f1 W ðtÞg þ
f1 f g: f1 xg b \

(18)

Here, the perturbation method is introduced to solve the non-linear equation. In Eq. (19), the small variant
½\1 kN\  can be regarded as the perturbation parameter term, because the variant
½\1 kN\  is small relative to ½\1 o2\ : The dynamic responses {1x} can be expanded in terms of a series of the small parameter
expressed as f1 xg ¼ f1 xð0Þ g þ
f1 xð1Þ g þ
2 f1 xð2Þ g þ ;

(19)

where superscripts denote the perturbation order. Then, the perturbed equations are evaluated as ð0Þ

f1 x€ g þ ½\1 o2\ f1 xð0Þ g ¼ f1 W ð0Þ g þ f1 f ð0Þ b g; ð1Þ

f1 x€ g þ ½\1 o2\ f1 xð1Þ g ¼ f1 f p ð1 xð0Þ Þg þ f1 f ð1Þ b g;

ð20Þ

ARTICLE IN PRESS B. Moon et al. / Mechanical Systems and Signal Processing 19 (2005) 1135–1151

1143

where {1fp} includes the non-linear stiffness term. {1W(0)} is the exciting Gaussian narrowband stationary random process. {1fb(0)} and {1fb(1)} are perturbed internal forces expressed by modal displacements. f1 f p ð1 xð0Þ Þg ¼ ½\1 kN\ f1 xð0Þ3 g: 2.4. Equations of motion of an assembled system The casing is modelled as the linear component ½2 Mf2 U€ c g þ ½2 Kf2 U c g ¼ f2 F E g þ f2 F b g;

(21)

where [2M] and [2K] are the mass and stiffness matrix, respectively. f2 F b g is the internal force vector. f2 F E gð¼ ½2 MfIgfX€ 0 gÞ is an earthquake external force term. After the eigenvalue analysis, Eq. (21) is changed into modal coordinate ðf2 U c g ½2 Ff2 xgÞ: The internal force is introduced in the equation because the linear component can be synthesised through the internal force with the other components. According to the previous section, Eq. (21) can be expressed with the approximated Gaussian stationary narrowband random process W(t). Even the casing component is a linear system. This component is perturbed in the same way as the non-linear component, because the higher order harmonic oscillation, which occurs in the non-linear component, is translated through the higher order perturbed equation as [3] ð0Þ

f2 x€ g þ ½\2 o2\ f2 xð0Þ g ¼ f2 W ð0Þ g þ f2 f ð0Þ b g; ð1Þ

f2 x€ g þ ½\2 o2\ f2 xð1Þ g ¼ f2 f ð0Þ b g;

ð22Þ

where {2fb(0)} and {2fb(1)} are the perturbed internal forces expressed by modal displacements. {2W(2)} is the external force term in the modal coordinate system. As an assembling component to apply the SSM, simple linear ball bearings are considered. Generally, there is a damping term in the bearing, but it is ignored in this study. The restoring force of the bearing is modelled as the linear term, where the force and displacement are expressed in matrix form as ½1 kb1 f1 U rb g ¼ f1 f b g;

½2 kb2 f2 U rb g ¼ f2 f b g;

(23)

where ½j kbj ðj ¼ 1; 2Þ are the terms of bearing coefficient. {1fb} and {2fb} are the internal force vectors of the non-linear component and linear component, respectively. {jUrb} is a relative displacement between the rotor and casing corresponding to the bearings. In order to solve the overall equation, the small parameter
is set equal to the perturbation parameter
of the nonlinear component. Then, the displacement is expressed as j ð1Þ fj U rb g ¼ fj U ð0Þ rb g þ
f U rb g

ðj ¼ 1; 2Þ:

(24)

Accordingly, by using Eq. (24), the internal force vectors can be perturbed as 1 ð0Þ f1 f b g ¼ f1 f ð0Þ b g þ
f f b g;

2 ð1Þ f2 f b g ¼ f2 f ð0Þ b g þ
f f b g:

(25)

ARTICLE IN PRESS B. Moon et al. / Mechanical Systems and Signal Processing 19 (2005) 1135–1151

1144

In order to synthesise the components, Eqs. (20), (22) and (25) are combined and rewritten in a matrix form. The equation of order
0 can be expressed as ð0Þ

ð0Þ

fx€ g þ ½K¯ fxð0Þ g ¼ fF ð0Þ ðtÞg;

(26)

where T 2 ð0Þ T 2 ð0Þ T T fxð0Þ g ¼ ff1 xð0Þ gT ; f1 U ð0Þ rb g ; f U rb g ; f x g g ; T 2 ð0Þ T 2 ð0Þ T T fF ðtÞð0Þ g ¼ ff1 W ð0Þ gT ; f1 f ð0Þ b g ;f fb g ;f W g g ð0Þ and ½K¯  is the stiffness matrix of the overall system, which is composed of all of the components. The equation of order
1 can be expressed as ð1Þ ð1Þ fx€ g þ ½K¯ fxð1Þ g ¼ fF ð1Þ ðxð0Þ Þg;

(27)

where T 2 ð1Þ T 2 ð1Þ T T fxð1Þ g ¼ ff1 xð1Þ gT ; f1 U ð1Þ rb g ; f U rb g ; f x g g ; T 2 ð1Þ T T T fF ð1Þ ðxð0Þ Þg ¼ ff1 f p ð1 xð0Þ ÞgT ; f1 f ð1Þ p g ; f f b g ; f0g g ; ð1Þ

where ½K¯  is the stiffness matrix of the overall system which is composed of all of the components. The overall system is obtained by assembling those component equations. The equation of order
0 is obtained as (

ð0Þ

f1 x€ g ð0Þ f2 x€ g

)

" þ

#(

½1 o2i  þ ½a1 

½a2 

½a3 

½2 o2i  þ ½a4 

fxð0Þ 1 g

(

) ¼

fxð0Þ 2 g

f1 f ð0Þ Z ðtÞg

)

f2 f ð0Þ Z ðtÞg

;

(28)

½a1  ¼ ½fb1 T ½1 kb1 ½fb1 ; ½a2  ¼ ½fb1 T ½2 kb1 ½fb2 ; ½a3  ¼ ½fb2 T ½1 kb2 ½fb1 ; ½a4  ¼ ½fb2 T ½2 kb2 ½fb2  2 ð0Þ where f1 f ð0Þ Z ðtÞg; f f Z ðtÞg are the external force terms. (

f1 f ð0Þ Z g

)

( ¼

f2 f ð0Þ Z g

½fa1 T f1 W ð0Þ gT ½fa2 T f2 W ð0Þ gT

) :

[fbi ] (i ¼ 1; 2) is the eigenvector matrix of the assembling region, which is derived from the eigenvector of each substructure corresponding to its bearing. The equation of order
1 is obtained as (

ð1Þ

f1 x€ g ð1Þ f2 x€ g

)

" þ

½1 o2i  þ ½a1 

½a2 

½a3 

½2 o2i  þ ½a4 

#(

fxð1Þ 1 g fxð1Þ 2 g

(

) ¼

f1 f ð1Þ Z g f2 f ð1Þ Z g

) ;

(29)

ARTICLE IN PRESS B. Moon et al. / Mechanical Systems and Signal Processing 19 (2005) 1135–1151

where (

f1 f ð1Þ Z g

)

( ¼

f2 f ð1Þ Z g

½fa1 T f½\1 kN\ f1 xð0Þ3 gg þ ½fb1 T f1 f ð1Þ b g

1145

)

½fa2 T f0g þ ½fb2 T f2 f ð1Þ b g

:

3. Stationary random properties of non-linear random vibration In this chapter, the response of random vibration of a non-linear system is solved statistically in an overall system. The earthquake is used as the excitation wave, which is regarded as the Gaussian stationary narrowband random process, by considering the strong motion duration of the earthquake. When the statistical properties of an earthquake are known, statistical properties of the system response can be obtained. After the eigenvalue analysis of the overall system with Eq. (28), the order
ð0Þ coordinate fZð0Þ g of the overall system is introduced for modal analysis as fxð0Þ g ½Ft fZð0Þ g;

(30)

where ½Ft  is the eigenvector matrix of the overall system. The damping of the system is assumed to be the proportional damping of the eigenvalue. Accordingly, the equation of the motion of order
ð0Þ is ð0Þ 2 ð0Þ _ ð0Þ Z€ ð0Þ i þ 2zoti Z i þ oti Zi ¼ W Zi

ði ¼ 1; 2; 3; . . . ; nÞ:

(31)

o2ti is the eigenvalue of the overall system. W ð0Þ Zi is the external force term in modal coordinates. By using linear random vibration theory, the solution of the linear differential equation (33) may be readily obtained. The response of
ð0Þ is Z 1 ð0Þ W ð0Þ (32) Zi ðtÞ ¼ Zi ðt  tÞhi ðtÞ dt; 0

where hi ðtÞ is the impulse response function of the linear system. Then, modal coordinates of the overall system fZð1Þ g for order
ð1Þ can be obtained as fxð1Þ g ½Ft fZð1Þ g:

(33)

By substituting Eq. (30) into Eq. (29), the equation of the motion of order
ð1Þ can be described as ð1Þ ð0Þ 2 ð1Þ _ ð1Þ Z€ ð1Þ i þ 2zoti Z i þ oti Zi ¼ f Zi ðZi Þ

ði ¼ 1; 2; 3; . . . ; nÞ;

(34)

ð0Þ 2 ð0Þ3 2 where f ð1Þ Zi ðZi Þð¼ bi Zi Þ is the external force term in modal coordinates. bi is the non-linear coefficient. Then, the response of
ð1Þ is Z 1 ð1Þ 2 Zð0Þ3 (35) Zi ðtÞ ¼ bi i ðt  tÞhi ðtÞ dt: 0

Accordingly, the response of a non-linear system can be evaluated by the perturbation theory as ð1Þ Zi ¼ Zð0Þ i þ
Zi :

(36)

ARTICLE IN PRESS 1146

B. Moon et al. / Mechanical Systems and Signal Processing 19 (2005) 1135–1151

ð1Þ The equations of Zð0Þ i and Zi can be used to compute various statistical properties of the response. Next, spectral analysis of the response Zi based on the series expansion is considered. For this, the covariance function is evaluated first. The covariance of the non-linear response, computed to the first order of
; can be obtained as ð0Þ ð0Þ ð1Þ ð1Þ ð0Þ EfZi ðt þ tÞZi ðtÞg ¼ E½Zð0Þ i ðt þ tÞZi ðtÞ þ
E½Zi ðt þ tÞZi ðtÞ þ
E½Zi ðt þ tÞZi ðtÞ:

(37)

This function can be expressed in other simple form with RZi Zi ðtÞ: RZiZi ðtÞ ¼ RZið0Þ Zið0Þ ðtÞ þ
RZið0Þ Zið1Þ ðtÞ þ
RZið1Þ Zið0Þ ðtÞ;

(38)

where RZið0Þ Zið0Þ ðtÞ is the covariance function of the linear part of the response, and RZið0Þ Zið1Þ ðtÞ and RZið1Þ Zið0Þ ðtÞ are the cross-covariance functions. After obtaining these terms, the covariance function of the non-linear response can be evaluated as  Z 1 2 3 2 2  1 ð0Þ ð0Þ S ðOÞ H i ðOÞ 
sZið0Þ dOSNi ðOÞ H i ðOÞ H i ðOÞ cos Ot ; (39) RZiZi ðtÞ ¼ 2 Ni 2 0 where H i ðOÞ is the conjugate function of H i ðOÞ: Then, the spectral density S ZiZi ðOÞ of the nonlinear response is obtained by taking the Fourier transform of the covariance function as 2 2 2 SZiZi ðOÞ ¼ Sð0Þ Ni ðOÞjH i ðOÞj ½1  6
bi sZið0Þ Re½H i ðOÞ;

(40)

where Re½H i ðOÞ is the real part of H i ðOÞ: The corresponding variance can be obtained from the covariance RZi Zi ðtÞ of the system by letting t ¼ 0; which is the same value of the mean-square response of the non-linear vibration. Z 1 2 2 2 fRZiZi ðtÞhi ðtÞg dt: (41) E½Zi ðtÞ ¼ RZiZii ð0Þ ¼ sZið0Þ ½1  6
bi 0

The spectral density is related to the stationary variance s2Zið0Þ ; which is the same as the meansquare value of the linear response. Examining s2Zi ; it appears that if the system is non-linear with light damping, weak non-linearity and the excitation random process is Gaussian stationary, then the response mean-square spectral density, correlation function, and mean-square value, can all be calculated from the knowledge of the mean-square spectral density SN ðOÞ of the excitation random process and the magnitude jHðOÞj of the frequency response. Table 1 Properties of the rotor system Rotor, casing length (mm) Rotor diameter (mm) Casing diameter (mm) Young’s modules of rotor and casing (N/m2) Density of rotor and casing (kg/m3) Bearing coefficient (N/m) Constrain coefficient (N/m)

800 16 50 2.1  1011 7.81  103 1.0  106 1.0  1010

ARTICLE IN PRESS B. Moon et al. / Mechanical Systems and Signal Processing 19 (2005) 1135–1151

1147

4. Numerical results Numerical calculation results of the non-linear stochastic system are obtained under the random excitation. Non-linear characteristics of the random response are observed by the proposed analytical method and those are compared with the results by the numerical simulation. 4.1. Model for analysis A non-linear rotor system, which is shown in Fig. 1, is considered. The rotor is considered as a uniform beam and the casing is also considered as a uniform beam approximately for the simplicity of calculation. The casing is constrained to a foundation at both ends of the casing. As a support, ball bearing is considered for the aircraft engine turbine or power plant turbine. The properties of the rotor system are tabulated in Table 1. The rotor is modelled by the 20 beam elements and the casing is also modelled by the twenty beam elements. The modal damping ratios of the system are given by z ¼ 0:05: To check the computing accuracy of SSM, the natural frequency of the rotor system is calculated by using FEM, and the results are compared with those calculated by using SSM. The rotor has 84 DOF because it has 21 modes and there are 4 variables per mode. The casing also has 84 DOF so that the total DOF are 168. When 40 modes are adopted in the rotor and casing, the error is less than 0.1%. To keep the accuracy of the lower natural frequencies, it is assumed that 20 modes of both components are sufficient for this analysis. The exciting force caused by the unbalance of the rotor is assumed as F u ðtÞ ¼ F 0 O2 cosðOt þ fÞ in x-direction, where F 0 and f are unbalance and phase angle, respectively. Rotating speed (O ¼ 540 rad=s) is assumed to be near first critical speed (ot1 ¼ 583 rad=s) of the rotor system. The unbalance of the rotor system is located at the middle of the shaft with a value of 0:044 kg m=ðrad=sÞ2 : 4.2. Linear response of the system Linear response of the system against the random excitation is examined. The responses are computed by SSM with the modal analysis procedure, which are shown in Fig. 4. The statistic properties of the system response are examined. As can be noticed from Fig. 4(a), the response includes the earthquake response of the system. Probability density of the system is shown in Fig. 4(b) in accordance with relative amplitude of the system. Average value of response is 0.000186 and standard deviation is 0.075741. These properties prove that the response is Gaussian stationary random process. In Fig. 4(c), the power spectrum of the response is shown. It can be observed that a typical earthquake response power spectrum is obtained, which has a low frequency component. The spectrum shows the fact that the response is in a short period and the predominant frequency is in a lower frequency where the peak frequency of the response is 28.27 rad/s. The simulated PSD shows a perfect with the analytical response of PSD S ZZ ðOÞ; which is obtained by
¼ 0:0: The spectral density is related to the stationary variance s2Z ð0Þ ¼ 0:0059; which is the same as the mean-square value of the linear response.

ARTICLE IN PRESS 1148

B. Moon et al. / Mechanical Systems and Signal Processing 19 (2005) 1135–1151

Fig. 4. Linear response of the system and its probabilistic properties.

4.3. Non-linear response of the system In this section, the statistical properties of the non-linear random vibration of the stochastic system verses the strong motion duration of excitation are investigated. Correlation, spectrum density and its variances, which are important to reliability analysis, are considered. Those properties of non-linear responses are calculated according to the procedure of the non-linear random vibration analysis, which is computed to the perturbation first order. Spectrum of harmonic excitation and earthquake excitation are obtained and applied to the system. If the dual inputs are uncorrelated, then the cross-spectral density function of earthquake and unbalance excitation terms in the equation of spectral density becomes zero. Then, the autocorrelation and the spectral density of the non-linear response are obtained theoretically. It is regarded that the non-linear response depends on the size of the perturbation parameter (e), which shows non-linear characteristics. To this end, two kinds of analysis are carried out, i.e. the perturbation parameters are set to
¼ 0:2 and 0.5.

ARTICLE IN PRESS B. Moon et al. / Mechanical Systems and Signal Processing 19 (2005) 1135–1151

1149

Fig. 5. Comparing the PSD of the system in accordance with various non-linear parameters. Table 2 Comparison of computing efficiency of the methods Analysis method

Variance/dev
¼ 0:0

Variance/dev
¼ 0:2

Variance/dev
¼ 0:5

Calculation time (s)

12 mode 20 mode 40 mode ELM DIM

0.00531/5.7 0.00522/3.9 0.00511/1.9 0.00520/3.5 0.00502

0.00472/4.6 0.00467/3.5 0.00461/2.2 0.00471/4.4 0.00451

0.00272/5.4 0.00269/4.2 0.00265/2.7 0.00267/3.4 0.00258

440 512 687 653 3370

In Fig. 5, the PSD of the non-linear responses of the system is shown with various non-linear parameters (=0.1, 0.4, 0.8). The response of SSM is calculated by taking 20 modes. Each PSD of the system shows a characteristic of the earthquake random process. Investigation of the PSD reveals that the PSD is smaller when the perturbation parameter becomes bigger, which is the non-linear response in terms of the non-linear characteristic. Fig. 5 shows another interesting phenomenon of the PSD function near the peak of the response and its shift from the frequency, known for non-linear systems. Variance of the non-linear response, which is an important value of the reliability analysis of the system, is evaluated from the non-linear response. Table 2 shows the variance of the non-linear response at the middle of the rotor by changing its numbers of adopting modes 12, 20 and 40. The values of the non-linear response are investigated according to the analytical methods for various values of non-linear parameter (
¼ 0:0; 0.2. 0.5). To prove the computing efficiency, those values are compared with the results of the direct integration method and Equivalent Linearised Method, which are obtained by same calculation conditions with the proposed method. The response of FEM is obtained by the direct integration of the equation of motion according to the Runge–Kutta method in which the equation of motion of the system is derived by the FEM, where the overall equation of motion of the system has 168 DOF.

ARTICLE IN PRESS 1150

B. Moon et al. / Mechanical Systems and Signal Processing 19 (2005) 1135–1151

Investigation of the variance reveals that the value shows a decrease with
in the spread of displacement about the equilibrium point when
¼ 0; as shown in Table 2. This is consistent with our intuition, which suggests that stiffer systems exhibit smaller displacements, and with the observation that the system stiffness increases with the non-linear parameter. This result also reveals that the variance displacement of a hardening spring-type non-linear system is always less than that for the corresponding linear system. Each analysis method stands for 12 modes, 20 modes and 40 modes (the response by the proposed method adopting 12 modes, 20 modes and 40 modes). ELM (the response by Equivalent Linearised Method), DIM (the response by direct integration method using the equation of FEM). Dev stands for the deviation of the calculation accuracy, which is noted as in Eq. (42) per percentage. Calculation accuracy of the proposed method is evaluated to show the effectiveness of the proposed method. The deviation of the calculation error of the analysis is defined as Deviation ¼

Value of DIM  Value of analytical methods  100 ð%Þ: Value of DIM

(42)

The deviations of the result to the direct integration method are 5.7%, 3.9% and 1.9% by the proposed method where the number of adopting modes is 12, 20 and 40, respectively, in the case of (
¼ 0:0). The deviations of the result to the direct integration method are 4.6%, 3.5% and 2.2% by the proposed method where the number of adopting modes is 12, 20 and 40, respectively, in the case of (
¼ 0:2). The deviations of the result to the direct integration method are 5.4%, 4.2% and 2.7% by the proposed method where the number of adopting modes is 12, 20 and 40, respectively, in the case of (
¼ 0:5). The deviations of the result to the direct integration method are 3.5%, 4.4% and 3.4% in the case of the Equivalent Linearised Method when (
¼ 0:0; 0.2, 0.5). The responses within 6% deviation error are obtained by the proposed analytical method. It is believed that the accuracy of the response within 6% deviation error for the vibration analysis of the system is the effective analytical method. Compared with the other conventional method, the proposed method is competitive in analysis accuracy by adopting only the low modes. Next, the calculation time is considered to verify the effectiveness of the proposed method. As a case, the calculation time for the variances of response in Table 2 is examined. The response of the direct integration method is obtained under the strong motion duration (=3–18 s) of excitation. The proposed method takes 440, 512 and 687 s to calculate the frequency response within 5% deviation error, where the number of adopting modes are 12, 20 and 40, respectively, when (
¼ 0:2), while the direct integration method takes 3370 s to compute the same response by using the personal computer Logix IBM Co. The Equivalent Linearised method takes 653 s to calculate the frequency response within 4.4% deviation error. As a result, it can be observed in this study that a drastic reduction in computational time can be obtained while retaining the accuracy of the solution. This is a critical factor in the analysis of the structural dynamics with a large number of DOF systems. It is believed that the proposed method can analyse the response of the complex system by keeping the accuracy of the solution compared with the other conventional method.

ARTICLE IN PRESS B. Moon et al. / Mechanical Systems and Signal Processing 19 (2005) 1135–1151

1151

5. Conclusions In this paper, the random vibration analysis method of a non-linear stochastic system was theoretically formulated using the perturbation method and the SSM. The actual random excitation is approximated to the corresponding Gaussian stationary random process for the analysis. The formulation is concerned with reducing the number of degrees of freedom for each component by modal substitution. All of the components are then assembled together and the random response of the overall system is analysed statistically against earthquake excitation near the 1st natural frequency of the system. It was shown that non-linear random responses could be efficiently calculated according to the selected number of vibration modes. Several statistical properties of the random responses that are of interest in non-linear random vibration applications are reviewed in accordance with the proposed method, which show the non-linear characteristics of the system. The variance value of the non-linear random response is obtained economically, which is important in evaluating the reliability of the system. The results reported herein provide a better understanding of the non-linear random vibration. Moreover, it is believed that those properties of the results can be utilised in the dynamic design of the non-linear stochastic system.

References [1] T. Iwatsubo, S. Kawamura, B. Moon, Non-linear vibration analysis of rotor system using substructure synthesis method (analysis with consideration of non-linearity of rotor), JSME International Journal Series III 41 (4) (1998) 727–733. [2] B. Moon, J. Kim, B. Yang, Non-linear vibration analysis of mechanical structure system using substructure synthesis method, KSME International Journal 13 (9) (1999) 620–629. [3] B. Moon, B. Kang, Non-linear vibration analysis of rotor system using substructure synthesis method (analysis with consideration of nonlinear sensitivity), JSME International Journal Series III 44 (1) (2001) 12–20. [4] A.H. Soni, V. Srinivasan, Seismic analysis of a gyroscopic mechanical system, Transactions of ASME, Journal of Vibration, Acoustics, Stress and Reliability in Design 105 (1983) 455–489. [5] V. Srinivasan, A.H. Soni, Seismic analysis of a rotor bearing system, Journal of Earthquake Engineering and Structure Dynamic 12 (1983) 287–311. [6] C.K. Rao, S. Mirza, Seismic analysis of high-speed rotating machinery, Nuclear Engineering and Design 111 (1989) 395–402. [7] M.F. Dimentgerg, M.N. Noori, Filtering property of a lightly damped system, subjected to transient excitation, Journal of Sound and Vibration 227 (2) (1999) 469–470. [8] M.F. Dimentgerg, Z. Hou, M.N. Noori, Spectral density of a non-linear single-degree-of-freedom system’s response to a white-noise random excitation: a unique case of an exact solution, International Journal of NonLinear Mechanics 30 (5) (1995) 673–676. [9] S. Bellizzi, R. Bouc, Spectral response of asymmetrical random oscillators, Probabilistic Engineering Mechanics 11 (1) (1996) 51–59. [10] S. Bellizzi, R. Bouc, Analysis of multi-degree of freedom strongly non-linear mechanical systems with random input: Part I: non-linear modes and stochastic averaging, Probabilistic Engineering Mechanics 14 (3) (1999) 229–244. [11] S. Bellizzi, R. Bouc, Analysis of multi-degree of freedom strongly non-linear mechanical systems with random input: Part II: equivalent linear system with random matrices and power spectral density matrix, Probabilistic Engineering Mechanics 14 (3) (1999) 245–256. [12] T.T. Soong, Mircea Grigoriu, Random vibration of mechanical and structural systems, pp. 64–65.