Application of fuzzy random finite element method on rotor dynamics
Qi Xu^{1} , Qian Zhao^{2} , Hongliang Yao^{3} , Bangchun Wen^{4}
^{1, 2, 3, 4}School of Mechanical Engineering and Automation, Northeastern University, Shenyang, China
^{3}Corresponding author
Journal of Vibroengineering, Vol. 16, Issue 4, 2014, p. 18541863.
Received 27 February 2014; received in revised form 13 April 2014; accepted 9 May 2014; published 30 June 2014
JVE Conferences
Fuzzy and stochastic characteristics of parameters exist widely in rotating machinery. To research the parameters characteristics is of great significance in rotor dynamics. Dynamic characteristics of rotor system are analyzed taking into account uncertain properties of fuzzy and stochastic coexisting. Fuzzy variables are transformed into stochastic variables based on information entropy theory. The Neumann stochastic finite element method based on Neumann expansion combined with Newmarkβ method is used in linear and nonlinear rotor system within the frame work of Monte Carlo simulation. Critical speed and dynamic response of fuzzy stochastic rotor systems are described by the proposed method. The results show that the Neumann stochastic finite element method has good applicability and efficiency in rotor dynamics.
Keywords: fuzzy, stochastic, rotor dynamics, Neumann expansion.
1. Introduction
Rotor dynamics is very important for responses prediction of rotating machineries. In the past 20 years, a lot of literatures were focused on numerical simulation of rotor systems, and dynamic models can be divided into three kinds: single disc Jeffcott rotor system with less degrees of freedom, multi discs and span rotor system with more degrees of freedom, complicated finite element model with nonlinear effects. And also, a lot of numerical methods have been applied, e.g., initial value methods such as RungeKutta method [1], Newmarkβ [2], Wilsonθ method; boundary value methods such as shooting method [3]. For example, Tejas H. Patel and Ashish K. Darpe [4] presented vibration response of a single disc rotor system with multi faults using RungeKutta method. Jerome Didier et al. [5] investigated nonlinear dynamic response of two discs rotor system with multi faults and uncertainties using finite element model. Nowadays finite element model is becoming popular increasingly to analyze complicated rotor systems [68], as it is capable to predict the static and dynamic behavior of rotor system based on its geometry and material characteristics, such as Mzaki Dakel et al. [9] modeled rotor system with hydrodynamic journal bearings using Timoshenko beam finite element model to study rotor nonlinear dynamics.
However, it is often difficult to define a reliable finite element model of rotor system with a number of uncertain physical properties. In fact, many components in rotor system are subject to uncertainty. There are two types of uncertainty in rotor system. One is probabilistic uncertainty. Uncertain parameters are described as random variables with known probability distributions. The other is fuzzy uncertainty; in contrast to probabilistic uncertainty, it is nonprobability as fuzzy data are not available.
Uncertainty has a great effect on dynamic behavior of rotor system greatly. Finite element model of a structure can be reliable by taking into account uncertainty. Many papers focused on the probabilistic uncertainty in rotor system, from random structures to random external forces [1012], such as Yanhong Ma et al. [13] considered support and connecting structure stiffness, phase and amount of rotor unbalance to present interval analysis method in rotor system. JeanJacques Sinou and Béatrice Faverjon [14] obtained dynamic response of a transverse crack in a rotating shaft by changing stiffness of the crack as random variables. Yazhao Qiu and Singiresu S. Rao [15] applied a fuzzy approach to analyze nonlinear rotorbearing system. Uncertainty analysis has been also developed quickly in other fields [1619], such as management science, engineering and technology, social economy, communication and transportation, finance and insurance. However, to the author’s knowledge, there is rarely work in rotor dynamics under the conditions of random and fuzzy coexisting.
In this paper, uncertain properties of fuzzy and random coexisting in rotor system are considered simultaneously to research rotor dynamics. Firstly, fuzzy variables are transformed into random variables based on information entropy theory, so Neumann expansion stochastic finite element method can be used in rotor dynamics. Secondly, Neumann expansion stochastic finite element method combined with Newmarkβ method is applied to analyze rotor system within the frame work of Monte Carlo simulation. Thirdly, critical speed and response are presented as examples in fuzzy random linear and nonlinear rotor systems by the proposed method. Efficiency of the method is verified by the examples.
2. Uncertainty in rotor system and entropy theory
2.1. Uncertainty in rotor system
Uncertainties exist widely in rotorbearing system due to randomness in material and geometric properties, or variant operation circumstance, and so on. External loads such as bearing reaction and rotating speed are all subject to variation; lubricant properties such as density and viscosity are varied by oil temperature; performance of components such as bearings and shafts is varied by wear and changes in operating conditions during their lifetimes.
Some of these uncertainties are statistical and probabilistic. For example, manufacturing and assembly tolerances exist in all mechanical parts and components, so dimensions of any mechanical part or component are probabilistic; eccentricity of wheel, clearance between the journal and the bearing, and so on are probabilistic uncertainty.
Some of the uncertainties do not have sufficiently reliable stochastic data, and they are associated with human error or limitation of professional knowledge. These uncertainties may be modeled by fuzziness. For example, sometimes it’s hard to determine whether boundary condition is simply supported or clamped supported, but can be described by using fuzzy terms, such as “nearly simple supported”, and “almost clamped supported”. In this case, it is fuzzy uncertainty.
2.2. Entropy theory and transformation rules
Information entropy is used to measure uncertainty of information, that is, entropy is a measure of uncertainty of a random variable [1113]. Probabilistic entropy for continuous random variable $X$ can be defined as:
where $H$ is entropy of $X$ and $p\left(x\right)$ is probability density function of $X$.
Fuzzy information can be measured by information entropy. It can be called fuzzy entropy. Nonprobabilistic entropy can be represented as:
where $\stackrel{\xb4}{G}$ is entropy of $X$ and $f\left(x\right)$ is membership function of $X$.
Fuzzy variable can be transformed into random variable based on entropy principle, that is, fuzzy entropy is equal to probabilistic entropy. The principle of transformation is that equivalent probabilistic entropy is equal to entropy of original fuzzy variable:
Normally, fuzzy variable is transformed into equivalent normal random variable. Assuming that normal random variable of mean is $m$ and standard deviation is $\sigma $. Probabilistic entropy of variable $X$ can be obtained by Eq. (1):
The equivalent standard deviation $\sigma $ can be given:
Fuzzy uncertainty can be transformed into any probabilistic distribution. Thus, equivalent normal random variable can be used to transform fuzzy structure to random structure for fuzzy stochastic finite element method applied in rotor dynamics.
Fig. 1. Three familiar types of fuzzy distributions: a) Triangular, b) Trapezoidal, c) $\Gamma $
a)
b)
c)
Three common types of fuzzy distributions are shown in Fig. 1. There are Triangular distribution, Trapezoidal distribution and $\mathrm{\Gamma}$ distribution, respectively. The expressions of the fuzzy distributions are as follows.
Triangular distribution:
Trapezoidal distribution:
$\mathrm{\Gamma}$ distribution:
$\stackrel{\xb4}{G}\left(x\right)=\left[\frac{k}{2\left({e}^{k{a}_{1}}+{e}^{k{a}_{1}}2\right)}\right]\left\{\mathrm{l}\mathrm{n}\left[\frac{k}{2\left({e}^{k{a}_{1}}1\right)}\right]1\right\}.$
3. Neumann expansion stochastic finite element method
Stochastic finite element method based on perturbation technique has been used widely to solve dynamic response of stochastic parameters structure with random excitation, but it’s not suitable for uncertain parameters structure which coefficient of variance are larger than 0.2.
The direct Monte Carlo simulation is suitable for problems with uncertain parameters which coefficient of variance are larger than 0.2, but the method is quite inefficient as large number of samples are required to guarantee accurate statistical result.
Neumann expansion stochastic finite element method, which based on the direct Monte Carlo method, can overcome the short coming of the direct Monte Carlo method, and is used to analyze uncertainty rotor system.
3.1. Neumann expansion theory and Newmarkβ method
3.1.1. Neumann expansion theory
Assuming ${\mathbf{K}}^{1}$ is the inverse matrix of $\mathbf{K}$, and matrix $\stackrel{}{\mathbf{K}}$ can be given $\stackrel{}{\mathbf{K}}=\mathbf{K}+\u2206\mathbf{K}$, so inverse matrix of $\stackrel{}{\mathbf{K}}$ is given by Neumann series:
where $\mathbf{I}$ is the identity matrix.
3.1.2. Newmarkβ method
Assuming $\mathbf{M}$, $\mathbf{C}$, $\mathbf{K}$ and $\mathbf{F}$ are mass, damping, stiffness and load matrix in rotor system, respectively. ${\mathbf{x}}_{0}$, ${\dot{\mathbf{x}}}_{0}$ and ${\ddot{\mathbf{x}}}_{0}$ are initial values. ${b}_{0}=1/\left(\beta {\u2206t}^{2}\right)$, ${b}_{1}=\delta /\left(\beta \u2206t\right)$, ${b}_{2}=1/\left(\beta \u2206t\right)\text{,}$${b}_{3}=1/2\beta 1\text{,}$${b}_{4}=\delta /\beta 1\text{,}$${b}_{5}=\u2206t/2\left(\delta /\beta 2\right)\text{,}$${b}_{6}=\u2206t\left(1\delta \right)\text{,}$ and ${b}_{7}=\delta \u2206t$ ($\delta \ge 1/2$, $\beta \ge 1/4{\left(1/2+\delta \right)}^{2}$) are the integral constants.
Effective stiffness matrix $\stackrel{~}{\mathbf{K}}$ is given $\stackrel{~}{\mathbf{K}}={b}_{0}\mathbf{M}+{b}_{1}\mathbf{C}+\mathbf{K}$.
Then triangular decomposition of $\stackrel{~}{\mathbf{K}}$ is represented $\stackrel{~}{\mathbf{K}}=\mathbf{L}\mathbf{D}{\mathbf{L}}^{T}$.
Effective load is given at $t+\u2206t$:
Then displacement vector is given:
Acceleration and velocity are calculated:
3.2. Neumann expansion Newmarkβ method
3.2.1. Linear system
Linear system dynamic equation is written when time is $t+\u2206t$:
The procedure using Newmarkβ method is as follows.
Firstly, effective stiffness matrix is given:
And then, for each time step, the effective load vector can be calculated:
Then triangular decomposition of the $\stackrel{~}{\mathit{K}}$ is represented $\stackrel{~}{\mathbf{K}}=\mathbf{L}\mathbf{D}{\mathbf{L}}^{T}$.
And then, node displacement vector at time $t+\u2206t$ can be solved:
Finally, the node acceleration and velocity can be calculated by:
The random stiffness matrix $\mathbf{K}$ can be decomposed into mean ${\mathbf{K}}_{0}$ and fluctuating deviator $\u2206\mathbf{K}$ by Neumann expansion theory, and $\mathbf{K}={\mathbf{K}}_{0}+\u2206\mathbf{K}$.
According to Neumann expansion method:
Only ${\mathbf{K}}_{0}$ is variable in $\mathbf{K}$ in linear system. Once inverse matrix of ${\mathbf{K}}_{0}$ is obtained, the Eq. (18) can be used iteratively to obtain inverse matrix of $\mathbf{K}$ without further decomposition and inverse, which can save a great deal of the CPU time.
Expansion series may be truncated when the series converges. In most cases, it can be expanded up to the third order to satisfy requirement of engineering.
3.2.2. Nonlinear system
Assuming dynamic equation of a weak nonlinear system (which is applied to most rotor systems) is given:
where $\mathbf{K}\left(\mathbf{x}\right)$ is nonlinear stiffness.
Effective stiffness matrix is given:
Assuming variable quantity of the stiffness matrix by Neumann expansion method is $\mathrm{\Delta}{\mathbf{K}}_{1}$, variable quantity of the stiffness matrix at the $i$th iteration in each time step is $\mathrm{\Delta}{\mathbf{K}}_{2}^{i}$:
For the first iteration:
For the $i$th iteration:
Calculation of the matrix decompose can be greatly reduced. $\mathbf{K}$ can be calculated by Neumann expansion theory.
4. Numerical examples
4.1. Example 1
A rotor system as shown in Fig. 2 has been considered for analysis. The system consists of two elastic supports, an elastic shaft and a rigid disc. There are 6 elements, 7 nodes and 28 degrees of freedom in finite element model of the rotor system. The supporting stiffness and damping ratio are 4.6×10^{7} N/m and 0.027, respectively. The rotating speed is 3000 rpm, and other data of the elements are shown in Table 1.
Fig. 2. Element model of rotor system
Dynamic equation of a $n$ nodes rotor system by Jalan A. K. and Mohanty A. R. [20] with finite element method can be written as:
where $\mathbf{M}$, $\mathbf{C}$, $\mathbf{G}$ and $\mathbf{K}$ are the mass, damping, gyroscopic moment and stiffness matrices of system, respectively, $\omega $ is angular velocity, $\mathbf{F}$ is exciting force vector.
Table 1. Data of elements
Element

1

2

3

4

5

6

Length / mm

80

80

10

10

80

80

Diameter / mm

10

10

80

80

10

10

First three order critical speeds of the rotor system are ${\omega}_{n1}=$65.68 Hz, ${\omega}_{n2}=$458.86 Hz, ${\omega}_{n3}=$1276.24 Hz, respectively.
Assuming the supporting stiffness obeys triangular distribution. The parameters are $a=$4.6×10^{7} N/m, ${a}_{1}=$4.0×10^{7} N/m and ${a}_{2}=$5.2×10^{7} N/m, respectively. Assuming the elastic modulus obeys normal distribution which mean value is 2.1×10^{11} Pa and variance is 0.05. Probability distributions of critical speed of the rotor system were carried out by Neumann expansion stochastic finite element method within 100000 times Monte Carlo simulation. Distributions of first three order critical speeds are shown in Fig. 3. Distribution parameters such as mean value and variance can be obtained by the distributions.
Fig. 3. Probability distribution of first three order critical speed in fuzzy and random rotor system: a) First order, b) Second order, c) Third order
a)
b)
c)
4.2. Example 2
Rotorstator rub fault often occurs on local position as small clearance between rotor and stator in rotor system. It can cause amplitude jumping phenomenon when rotor system runs. It is a distinct nonlinear phenomenon, and can be explained by nonlinear stiffness in rotor system.
Fig. 4. Element model of nonlinear rotor system
Assuming the nonlinear stiffness are ${\stackrel{\xb4}{k}}_{x}={\stackrel{\xb4}{k}}_{y}=$5×10^{7} N/m at node 3 in rotor system as shown in Fig. 4. Nonlinear force can be given:
Amplitudefrequency response curve was carried out by Newmarkβ method as shown in Fig. 5. It shows that the system occurred amplitude jumping phenomenon. The jumping point is $\mathrm{a}$ when the rotor system is runup, and the jumping point is $b$ when the rotor system is rundown.
Fig. 5. Amplitude jumping phenomenon in nonlinear rotor system
Fig. 6. Probability distributions of jumping point
Assuming the supporting stiffness obeys triangular distribution. The parameters are $a=$5×10^{7} N/m, ${a}_{1}=$2×10^{7} N/m and ${a}_{2}=$8×10^{7} N/m, respectively. Assuming eccentric mass obeys normal distribution which mean value is 1×10^{2} kg/m and variance is 0.05. Probability distribution of the jumping point was calculated by Neumann expansion stochastic finite element method within 1000 times Monte Carlo simulation as shown in Fig. 6. Distribution parameters such as mean value and variance can be obtained by the distributions.
4.3. Example 3
Rub fault often occurs in compressor rotor system. Finite model of a compressor medium pressure cylinder rotor system is shown in Fig. 7. There are 23 elements, 24 nodes and 96 degrees of freedom in the finite element model. The supporting stiffness and damping are 12×10^{8} N/m and 11×10^{5} N$\bullet $s/m, respectively. Local rub fault occurs at node 12, and mass eccentricity is at nodes 9 and 16.
Fig. 7. Element model of compressor medium pressure cylinder rotor system
First three order critical speeds of the rotor system are ${\omega}_{n1}=$58.42 Hz, ${\omega}_{n2}=$182.26 Hz, ${\omega}_{n3}=$72.98 Hz, respectively. Amplitudefrequency response curve was carried out by Newmarkβ method as shown in Fig. 8.
Fig. 8. Amplitudefrequency response curve of the system
Assuming the supporting stiffness obeys triangular distribution. The parameters are $a=\text{1}\text{2}\text{\xd7}\text{10}\text{8}\text{}\text{N/m,}$${a}_{1}=\text{8}\text{\xd7}\text{10}\text{8}\text{}\text{N/m}$ and ${a}_{2}=\text{1}\text{6}\text{\xd7}\text{10}\text{8}\text{}\text{N/m,}$ respectively. Assuming elastic modulus obeys normal distribution which mean value is 2.1×10^{11} Pa and variance is 0.05. Probability distributions of critical speed, peek and peek frequency of the rotor system were carried out by Neumann expansion stochastic finite element method. 100000, 10000 and 1000 times Monte Carlo simulation were calculated, respectively. Distributions of first three order critical speeds, peek and peek frequency are shown in Figs. 9 and 10. Distribution parameters such as mean value and variance can be obtained by the distributions.
Fig. 9. Probability distribution of first three order critical speed in compressor medium pressure cylinder rotor system: a) First order, b) Second order, c) Third order
a)
b)
c)
Fig. 10. Probability distributions of peek and peek frequency: a) peek, b) peek frequency
a)
b)
5. Conclusions
The Neumann stochastic finite element method based on Neumann expansion combined with Newmarkβ method is used in rotor dynamic with fuzzy and random coexisting. The analysis indicates that dynamic characteristics such as critical speed, amplitude jumping phenomenon and peek frequency are affected in uncertainty system. The proposed method can be used in uncertainty rotor system of linear or nonlinear with the frame work of Monte Carlo simulation, and overcome the limit that coefficient of variance cannot be larger than 0.2. Meanwhile, a large of computational complexity and time of matrix decomposition and inverse can be reduced by the proposed method. The examples show that the method applied in rotor dynamic is effective.
Acknowledgements
The authors would like to gratefully acknowledge the National Basic Research Program of China (2011CB706504), the National Natural Science Foundation of China (51005042) and the Fundamental Research Funds for the Central Universities of China (N120403007).
References
 ChangJian C. W., Chen C. K. Nonlinear dynamic analysis of rubimpact rotor supported by turbulent journal bearings with nonlinear suspension. International Journal of Mechanical Sciences, Vol. 50, Issue 6, 2008, p. 10901113. [Search CrossRef]
 Lee A. S., Kim B. O., Kim Y. C. A finite element transient response analysis method of a rotorbearing system to base shock excitations using the statespace Newmark scheme and comparisons with experiments. Journal of Sound and Vibration, Vol. 297, Issue 35, 2006, p. 595615. [Search CrossRef]
 Musznyska A. Rotor Dynamics. Taylor & Francis, New York, 2005. [Search CrossRef]
 Patel T. H., Darpe A. K. Vibration response of a cracked rotor in presence of rotorstator rub. Journal of Sound and Vibration, Vol. 317, Issue 35, 2008, p. 841865. [Search CrossRef]
 Didier J., Sinou J. J., Faverjon B. Study of the nonlinear dynamic response of a rotor system with faults and uncertainties. Journal of Sound and Vibration, Vol. 331, Issue 3, 2012, p. 671703. [Search CrossRef]
 Xiang J. W., Chen D. D., Chen X. F., He Z. J. A novel waveletbased finite element method for the analysis of rotorbearing systems. Finite Element in Analysis Design, Vol. 45, Issue 12, 2009, p. 908916. [Search CrossRef]
 Wang Z. L., Yu X. L., Liu F. L., Feng Q. K., Tan Q. Dynamic analyses for the rotorjournal bearing system of a variable speed rotary compressor. International Journal of Refrigeration, Vol. 36, Issue 7, 2013, p. 19381950. [Search CrossRef]
 Taplak H., Parlak M. Evaluation of gas turbine rotor dynamic analysis using the ﬁnite element method. Measurement, Vol. 45, Issue 5, 2012, p. 10891097. [Search CrossRef]
 Dakel M., Baguet S., Dufour R. Nonlinear dynamics of a supportexcited flexible rotor with hydrodynamic journal bearings. Journal of Sound and Vibration, Vol. 333, Issue 10, 2014, p. 27742799. [Search CrossRef]
 Young T. H., Shiau T. N., Kuo Z. H. Dynamic stability of rotorbearing systems subjected to random axial forces. Journal of Sound and Vibration, Vol. 305, Issue 3, 2007, p. 467480. [Search CrossRef]
 Leng X. L., Meng G. A., Zhang T., Fang T. Bifurcation and chaos response of a cracked rotor with random disturbance. Journal of Sound and Vibration, Vol. 299, Issue 3, 2007, p. 621632. [Search CrossRef]
 Qiu Y. Z., Rao S. S. A fuzzy approach for the analysis of unbalanced nonlinear rotor systems. Journal of Sound and Vibration, Vol. 284, Issue 12, 2005, p. 299323. [Search CrossRef]
 Ma Y. H., Liang Z. C., Chen M., Hong J. Interval analysis of rotor dynamic response with uncertain parameters. Journal of Sound and Vibration, Vol. 332, Issue 16, 2013, p. 38693880. [Search CrossRef]
 Sinou J. J., Faverjon B. The vibration signature of chordal cracks in a rotor system including uncertainties. Journal of Sound and Vibration, Vol. 331, Issue 1, 2012, p. 138154. [Search CrossRef]
 Qiu Y. Z., Rao S. S. A fuzzy approach for the analysis of unbalanced nonlinear rotor systems. Journal of Sound and Vibration, Vol. 284, Issue 12, 2005, p. 299323. [Search CrossRef]
 Trawińskia K., Cordóna O., Quirinc A., Sánchezd L. Multiobjective genetic classifier selection for random oracles fuzzy rulebased classifier ensembles: How beneficial is the additional diversity? KnowledgeBased Systems, Vol. 54, 2013, p. 321. [Search CrossRef]
 Shapiro A. F. Modeling future lifetime as a fuzzy random variable. Insurance: Mathematics and Economics, Vol. 53, Issue 3, 2013, p. 864870. [Search CrossRef]
 Parka J. P., Jeong J. U. On random fuzzy functional differential equations. Fuzzy Sets and Systems, Vol. 223, 2013, p. 8999. [Search CrossRef]
 Wang W. J., Liu D., Liu X., Pan L. Fuzzy overlapping community detection based on local random walk and multidimensional scaling. Physica A: Statistical Mechanics and its Applications, Vol. 392, Issue 24, 2013, p. 65786586. [Search CrossRef]
 Jalan A. K., Mohanty A. R. Model based fault diagnosis of a rotorbearing system for misalignment and unbalance under steadystate condition. Journal of Sound and Vibration, Vol. 327, Issue 35, 2009, p. 604622. [Search CrossRef]