Uncertainty representation and quantification for a nonlinear rotor/stator system with mixed uncertainties

Lechang Yang1 , Jianguo Zhang2 , Yanling Guo3

1, 2School of Reliability and Systems Engineering, Beihang University, Beijing 100191, China

1, 2Science and Technology on Reliability and Environmental Engineering Laboratory, Beihang University, Beijing 100191, China

3School of Automation Science and Electrical Engineering, Beihang University, Beijing 100191, China

1Corresponding author

Journal of Vibroengineering, Vol. 18, Issue 7, 2016, p. 4836-4851. https://doi.org/10.21595/jve.2016.17340
Received 28 June 2016; received in revised form 4 September 2016; accepted 25 October 2016; published 15 November 2016

Copyright © 2016 JVE International Ltd. This is an open access article distributed under the Creative Commons Attribution License, which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.
Creative Commons License

A rotor-to-stator coupled system usually exhibits complicated dynamic behaviors due to its nonlinear nature. Moreover, the inherent uncertainty (aleatory uncertainty) and many undetermined factors either brought by manufacturing process or due to the lack of knowledge (epistemic uncertainty) make the analysis of system response a challenging task. Existing studies on rotor uncertainties are mostly focused on the stochastic variables, yet pay less attention to other forms of uncertain variables such as intervals. However, some physical parameters (e.g. friction coefficient) can be hardly assigned one specific probability distribution and often available in interval forms. To deal with this, the concept of likelihood is extended from classical discrete point value to interval variable in the presence of mixed uncertainties. A likelihood-based approach is carried out for the mixed uncertainties representation and quantification. In addition, a new single loop sampling algorithm is developed to reduce the computation cost. This framework could be applied in the field of industry manufacturing and mounting, especially take effect in risk assessment and product maintaining. A series of numerical cases are demonstrated for validation and comparison.

Keywords: nonlinear rotor/stator system, mixed uncertainties, likelihood-based representation, uncertainty quantification.

1. Introduction

Motivated by the requirement to enhance the efficiency in engineering, the speed of rotating machinery is continuously increasing while the gap between the rotor and stator is significantly reduced. Affected by various uncertain factors, the rotor/stator rubbing is not rare in practical engineering. Unexpected rub induced by instability may lead to catastrophic consequence.

In the deterministic research field, a large amount of work has been done to better understand the rub-related phenomenon and get deep insights into its physical nature. Black [1] first discovered the existence condition of the synchronous full annular rub solution in early 1680s. Later Muszynska [2, 3] reported many possible responses of rotor/stator system. Through the numerical and experimental investigations, many typical contact responses as well as nonlinear dynamic behaviors have been captured, for instance, the jump phenomenon and the synchronous full annular rubs [4, 5], the partial rubs in sub- and super-synchronous whirl [6-8], the partial rubs in quasi-periodic whirl [9, 10], the chaotic motion [11-13] as well as dry whip [14-17]. In the last decade, Jiang etc. [18, 19] discussed the physical mechanism of dry whip and analytically derived the onset/existence conditions for a modified Jeffcott rotor.

However, the parameter uncertainty is inherent in most practical engineering and its effect on system response should not be overlooked. This problem has been realized in recent years and addressed in a series of publications [20-23]. These research works have made great contributions to the field of design and assessment of rotor/stator system considering parameter uncertainty. However, most of these monographs are concentrating on the parameters stochastic nature. From the perspective of uncertainty theory, the uncertainty that they are dealing with is pure aleatory uncertainty. The global response of a rotor/stator system under mixed uncertainty is still unclear. This issue is challenging in dynamic analysis or reliability evaluation and has not been sufficiently addressed yet. To fill the research gap and establish a more generic framework, this article proposes a likelihood-based approach for uncertain dynamic analysis or reliability estimation. Limitations embedded in traditional methods have been eliminated by extenging the concept of likelihood from traditional sparse point data to interval information. A single loop sampling algorithm is developed to simplify the computing process.

This paper is organized as follows, the modified Jeffcott rotor/stator system and its corresponding governing equations are presented in section 2, along with the analyses of dynamics and stability under ideal circumstance (deterministic) in section 3. The basic principle and detailed algorithm of the proposed likelihood-based approach are illustrated in section 4. Several numerical cases are demonstrated for validation and comparison purpose in section 5. The benefits and possible experimental verification are discussed in section 6. Finally, some concluding remarks are drawn.

2. Mathematical model and motion equation

A modified Jeffcott rotor/stator system is taken into study shown as Fig. 1, which is comprised of an eccentric rotor and a rigidly fixed stator with elastic surface. The rotor consists of a massless shaft carrying a disk with mass m mounted at the middle of the span. The disk mass center deflects its geometrical center with a distance e, namely, the rotor has an eccentricity e. This leads a mass imbalance when the rotor rotates at a speed ω. The stiffness of rotor shaft is denoted as kr and the stator is approximated modeled by radial springs with stiffness ks. There is an initial clearance r0 between the rotor and stator.

Fig. 1. a) Schematic diagram of the rotor/stator system, b) applied forces of rotor during whirling

 a) Schematic diagram of the rotor/stator system, b) applied forces of rotor during whirling


 a) Schematic diagram of the rotor/stator system, b) applied forces of rotor during whirling


The motion equation is given as:

m x ¨ + c x ˙ + k r x + s i g n k s 1 - r 0 r x - f c y = m e ω 2 c o s ω t , m y ¨ + c y ˙ + k r y + s i g n k s 1 - r 0 r f c x + y = m e ω 2 s i n ω t ,
s i g n = 0 ,         x 2 + y 2 < r 0 , 1 ,         x 2 + y 2 > r 0 ,

where the damping ratio, Coulomb friction coefficient are denoted by c and fc respectively. Eq. (1) can be rewritten in the non-dimensional form and formulated as:

X ¨ + 2 ζ X ˙ + β X + s i g n 1 - R 0 R X - f c Y = Ω 2 c o s Ω τ , Y ¨ + 2 ζ Y ¨ + β Y + s i g n 1 - R 0 R f c X + Y = Ω 2 s i n Ω τ ,         s i g n = 0 ,         X 2 + Y 2 < R 0 , 1 ,         X 2 + Y 2 > R 0 ,

where the non-dimensional variables are defined as:

X = x e ,           Y = y e ,           R 0 = r 0 e ,           R = r e ,           ω 2 2 = k s m ,
c m = 2 ζ ω 2 ,           Ω = ω ω 2 ,           β = k r k s ,           τ = ω 2 t .

It is apparent that the motion equation is linear when the rotor is separate from the stator (i.e. sign= 0). However, if the displacement of the rotor center exceeds initial clearance, the rotor would get in contact with the stator and the motion equation becomes nonlinear. This rotor/stator system has complicated dynamic behaviors with nonlinear characteristics. The various responses of the rotor/stator system (periodic full annular rub, quasi-periodic partial rub and dry whip, see in Fig. 2) with different given parameters will be discussed in the next section.

Fig. 2. Rotor orbit for different parameters: a) full annular rub, b) partial rub, c) dry whip

Rotor orbit for different parameters: a) full annular rub, b) partial rub, c) dry whip


Rotor orbit for different parameters: a) full annular rub, b) partial rub, c) dry whip


Rotor orbit for different parameters: a) full annular rub, b) partial rub, c) dry whip


3. Deterministic system response and stability analysis

Since our goal is to better understand the responses of rotor/stator system considering uncertainty, it is necessary to study the dynamics behaviors under deterministic condition. There are multi-types of responses of the rotor/stator system which are governed by different motion equations depending on the input parameters. For instance, if the deflection of the rotor center is less than the gap between rotor and stator, the rotor is not in contact with stator and the motion equation becomes linear. When the rotor gets full contact with the stator, it is the full annular rub that has been observed in the experiment [17, 24, 25]. However, if the rotation speed continuous goes up, the rotor may lose its stability, the system response becomes quasi-periodic and the motion is partial rub. In other cases, the friction direction is changed due to the high rotation speed or crude contact surface, which leads to a backward whirl motion (dry whip). Relative research works could refer [14-16]. For mathematical convenience, the motion equation is rewritten in the vector form as:

d W = L W + s i g n N L ,


W = X Y X ' Y ' ,           L = 0 0 1 0 0 0 0 1 - β 0 - 2 ζ 0 0 - β 0 - 2 ζ ,           N L = 0 0 1 - R 0 / R X - f c Y + Ω 2 c o s Ω τ 1 - R 0 / R f c X + Y + Ω 2 s i n Ω τ .

The steady-state periodic solution is in the form as X0=RcosΩτ+ψ, Y0=RsinΩτ+ψ. Eq. (3) is linearized about the steady-state solution (X0, Y0) as:

d δ W _ = J δ W _ ,

where d is a derivative operator and J is the Jacobian matrix. If we use W to represent a perturbation of the steady solution (X0, Y0), the stability of W is able to reflect the stability of solution (X0, Y0). Hence, the stability of the corresponding solution of Eq. (4) could be utilized to determine stability of solutions in Eq. (3) and thus, the stability for both linear solution and nonlinear solution [5].

3.1. Deterministic boundary for linear solution

For the linear scenario, the rotor will not contact with the stator, namely, sign= 0. Thus, the govern equation reduce to:

d W = L W .

The steady-state periodic solution is presented in the form as X=RcosΩτ+ψY=RsinΩτ+ψ. The amplitude and phase angle are given by:

R = Ω 2 β - Ω 2 2 + 4 ζ 2 Ω 2 ,           t a n ψ = 2 ζ Ω β - Ω 2 .

However, this solution is only eligible when the amplitude R is less than clearance R0. Let RR0, the parameter equation for no rub periodic solution is derived as Eq. (7):

R 0 2 - 1 Ω 4 + R 0 2 4 ζ 2 - 2 β Ω 2 - 4 R 0 2 ζ Ω + R 0 2 β 2 0 .

On the other hand, the associated characteristic equation for Jacobian matrix J is expressed as:

λ 4 + 4 ζ λ 3 + 4 ζ 2 + 2 β λ 2 + 4 ζ β λ + β 2 = 0 .

By adopting the Routh-Hurwitz condition, it is found that the solution is always stable since 4ζ2β>0. Thus, the two real roots Ωl and Ωu, which solved from Eq. (7), can be used to determine the stable boundary of no-rub motion.

3.2. Deterministic boundary for nonlinear solution

When the rotor get in contact with the stator, i.e. sign= 1 the governing equation becomes nonlinear. For this scenario, the Jacobian matrix J is time-dependent and is in following form:

J = 0 0 1 0 0 0 0 1 - 1 - β + R 0 R 1 - c o s 2 φ + f c R 0 R s i n φ c o s φ f c 1 - R 0 R + s i n 2 φ - R 0 R s i n φ c o s φ - 2 ξ 0 - f c 1 - R 0 R + c o s 2 φ + R 0 R s i n φ c o s φ - 1 - β + R 0 R 1 - s i n 2 φ - f c R 0 R s i n φ c o s φ 0 - 2 ξ .

To study the stability of one specific solution, corresponding transformation is taken as:

δ W = T δ U ,


T = c o s φ - s i n φ s i n φ c o s φ c o s φ - s i n φ s i n φ c o s φ ,

thus, Eq. (11) is derived by substituting Eq. (10) into Eq. (4):

d T δ U = J T δ U ,
d T δ U + T d δ U = J T δ U ,
δ U = T - 1 J T - d T δ U .

Define J_=T-1JT-dT and the J is solved as:

J _ = 0 Ω 1 0 - Ω 0 0 1 - 1 - β f c 1 - R 0 / R - 2 ζ Ω - f c - 1 - β + R 0 / R - Ω - 2 ζ .

Now, the time-dependent component T is separated from the Jacobian matrix J and the time-independent matrix J can be used to evaluate the stability of solution. The characteristic equation is written as:

λ 4 + b 3 λ 3 + b 2 λ 2 + b 1 λ + b 0 = 0 ,

where λ denotes the eigenvalues of J and coefficients are determined as:

b 3 = 4 ζ ,           b 2 = 4 ζ 2 + 2 Ω 2 + 2 1 + β - R 0 / R ,
b 1 = 4 f c Ω + ζ + β ζ + Ω 2 ζ - 2 f c Ω + ζ R 0 / R ,
b 0 = Ω 2 Ω 2 + 4 ζ 2 - 2 + 1 + β 2 + f c 2 + Ω 2 - f c 2 - 1 - β R 0 / R .

As these eigenvalues are related to the amplitude R, and R should be first solved and substituted into Eq. (13). Then the Routh-Hurwitz criteria can be adopted to evaluate the stability of nonlinear solution. For more details, please refer [5]. However, since our goal is not only to judge the stability of one given solution but also to determine the safe region in the parameter space, the parameter boundary is analyzed from the bifurcation perspective.

3.2.1. Saddle-node bifurcation

Based on the saddle-node bifurcation theory, there should be one zero eigenvalue of the Jacobian matrix. By substituting λ= 0 into characteristic equation, it leads to:

Ω 2 Ω 2 + 4 ζ 2 - 2 + 1 + β 2 + f c 2 + Ω 2 - f c 2 - 1 - β R 0 / R = 0 .

By eliminating the amplitude R through symbolic calculation, a twelve-order polynomial of rotation speed Ω is derived. With the aid of numerical simulation, saddle-node bifurcation boundaries can be solved and depicted on the parameter region.

3.2.2. Hopf bifurcation

According to the Hopf bifurcation theory, there should be one pair of conjugate purely imaginary eigenvalues (e.g. λ=iω or λ=-iω) for the Jacobian matrix. By substituting λ=iω (or λ=-iω) into the characteristic equation, Eq. (15) is derived as:

ω 4 - b 2 ω 2 + b 0 + - b 3 ω 3 + b 1 ω i = 0 .

Eliminating ω, it leads to:

b 1 2 - b 1 b 2 b 3 + b 0 b 3 2 = 0 , b 1 b 3 > 0 .

By Substituting b0-b3 in Eq. (13) into Eq. (16) and after manipulation, Eq. (17) is derived as:

c 2 R 2 + c 1 R + c 0 = 0 ,

where the coefficient c0-c2 is given as:

c 2 = 16 4 ζ 2 1 + β - f c 2 ζ 2 + Ω 2 ,           c 1 = 16 R 0 f c β - 2 ζ 2 Ω 2 + ζ 2 ,

After eliminating R, a twelve-order polynomial of rotation speed Ω is also obtained and the boundary of Hopf bifurcation can be numerically solved.

4. Likelihood-based representation of epistemic uncertainty due to sparse point data and interval data

In previous section, the system response in ideal circumstance, i.e. the deterministic boundary in parameter region, has been studied. In this section, multi-source uncertain factors within the system are taken into account and a likelihood-based approach dealing with mixed uncertainty is developed.

4.1. Description of uncertainty

In practical engineering, many uncertain factors (e.g. physical uncertainty, model uncertainty) are induced either by biased experiential judgments or by lack of knowledge. To quantify them, an unified measurement need to be built. According to the recently developed uncertainty theory, the uncertainty can be mainly classified into two groups, i.e. aleatory uncertainty and epistemic uncertainty. This distinction for uncertainty has been recognized by the NNSA and emphasized in the NAS/NRC report about QMU (Quantification of Margins and Uncertainties). For the basis and rationale of uncertainty theorem, please refer relevant monographs or publications [26-32].

The aleatory uncertainty arises from an inherent randomness in the properties or behavior of the system under study. On the other hand, the epistemic uncertainty derives from a lack of knowledge about the appropriate value to use for a quantity that is assumed to have a fixed value in the context of a particular analysis [26]. The aleatory uncertainty is supposed to be irreducible and usually represented in the stochastic form. By contrast, the epistemic uncertainty can be reduced when new information becomes available, it is represented in more complicated forms as both probabilistic and non-probability variable may coexist. (e.g. interval or point data). In this paper, both aleatory and epistemic uncertainties are taken into consideration.

4.2. Construction to the likelihood of uncertain interval data

In this subsection, the concept of likelihood function is first reviewed and extended from classical sparse point data (deterministic) to interval information (uncertain), which serves the foundation of the likelihood-based approach.

Consider a stochastic variable X follows one specific distribution (e.g. Normal) with a series of collected point data, its PDF (Probability Density Function) is denoted as px|Θ where Θ corresponds to the parameter set. The likelihood function LΘ|x is defined as the probability of a series of given data set conditioned on distribution parameter Θ [30, 33, 34]. Theoretically speaking, the probability value for any single discrete point is zero for a continuous density function, the likelihood function is obtained by considering the probability of an infinitesimally small interval around xi based on mean value theorem as Eq. (18) for practical use in previous work [30, 33, 34]:

L Θ P X x i - ε 2 , x i + ε 2 | Θ = x i - ε / 2 x i + ε / 2 f X x | Θ d x = ε f X x i | Θ f X x i | Θ .

The parameter Θ can be evaluated by maximizing the right hand side of Eq. (18). This estimation approach is popularly known and widely adopted as the maximum likelihood estimate (MLE).

Similar principle can be extended to bounded interval data, thus the expression for likelihood of the parameters Θ for interval [a, b] is:

L Θ P a x b | Θ = a b f X x | Θ d x = F X b | Θ - F X a | Θ ,

where the FX  corresponds to the CDF (Cumulative Density Function) of variable X. For multiple input interval information, the combined likelihood function is expressed as:

L Θ i = 1 n F X b i | Θ - F X a i | Θ .

It is the CDF rather than PDF that will be used in the likelihood calculation for interval data, compared with that for point data. This is the main difference between the construction of likelihood for conventional point data and interval data.

4.3. Representation of epistemic uncertainty with additional probabilistic information

Sometimes, additional probabilistic information is available for the uncertainty variable besides the given intervals and/or sparse point data. For instance, the rotor eccentricity e may be given by a series of intervals based on expert judgments and some sparse point collected from the measurements of similar product in this case. On the other hand, from the physical background perspective, the uncertainty of eccentricity e mainly arise from the error/deviation accumulated during the procedure of manufacturing, mounting, etc., which generally follows the normal distribution (statistics parameter is unknown).

For this scenario, if variable X is supposed to follow one specific distribution, the combined likelihood function can be expressed as Eq. (21). Since the comprehensive likelihood is expressed as Eq. (21), the distribution of parameter Θ can be derived utilizing the MLE. However, this calculation may be cumbersome because the cumulative density function of variable X is not always in analytical form (e.g. the CDF of normal distribution). In this paper, instead of computing statistics parameter Θ by maximizing the likelihood, a Bayesian-based approach is adopted to estimate the full likelihood. The idea that the entire likelihood function, rather than merely its maximizer, should be used for inference was emphasized by Barnard et al. [35]:

L Θ i = 1 n f X x i | | Θ i = 1 n F X b i | Θ - F X a i | Θ .

Assigning an appropriate prior distribution, the parameter Θ is estimated as Eq. (22) by adopting Bayesian theory:

p Θ | D = l Θ | D π Θ l Θ | D π Θ d Θ ,

where πΘ denote the joint prior distribution of parameter vector Θ, while pΘ|D is the estimated probability distribution of Θ. lΘ|D is the combined likelihood function which has the same form as Eq. (21).

It should be noted that the parameter Θ estimated in Bayesian method is a probability distribution, compared with a deterministic value by MLE method. Hence, a family of conditioned probability distribution for X is derived due to the uncertainty within the statistics parameter θ. A double loop Monte Carlo (or called as second-order Monte Carlo) method (Table 1 shows) is usually adopted for analysis, but this method may not be affordable for some complicated cases. To illustrate, suppose M samples of θ are generated in the outer loop, each of which corresponds to a distribution of pθ, and N samples of X are drawn for each sample θ in the inner loop, the computation cost is total M×N. To get a satisfactory result, 104 samples is usually required for M and N, thus the number of total samples is incredible large.

Table 1. Double loop Monte Carlo method

Suppose the model output Y is given by a deterministic function of Y=fX
To get an estimated result of Y:
Outer loop:
M samples of parameter θ are generated from pθ
Inner loop:
For a given θ
N samples of X generated from pXθ
End loop
Output Y is calculated as Y=fX for a family of distributions of variable X
End loop
The PDF of model output Y is estimated based on the M×N samples (such as kernel density estimate)

Considering this, we discard the double loop sampling approach and construct the PDF of variable X as Eq. (23) to replace the family of distribution for pX|θ with a unique distribution pX:

p x = p x | Θ p Θ d Θ .

The unconditional PDF pX can be calculated by integrating the conditional PDF px|Θ over the entire parameter space of Θ.

4.4. Representation of epistemic uncertainty without additional probabilistic information

In subsection 4.3, a methodology of representing the epistemic uncertainty with additional probability information is present. In other cases, additional probabilistic information is unobtainable, only single/multiple interval(s) is available. For instance, it is difficult to determine the probability distribution for the friction coefficient f.

Hence, a methodology to deal with pure interval uncertainty is required. To begin with, a sole interval [a, b] is assigned to the uncertain variable X, i.e. Xa,b. From the probabilistic perspective, it can be interpreted as X~Ua,b since there is an equal chance for variable X taking any value within the range of [a, b]. It can also be interpreted as an optimization problem as following

Object: fx

Subject: (1) fx0, (2) abfrxdx=1, (3) fxp.

The solved PDF fx is the corresponding probabilistic description for the given interval [a, b].

For multiple N input intervals xai,bi, 1iN, there are total 2N-1 subintervals dj,dj+1, amindjdj+1bmax, where amin is the lower bound and bmax is the upper bound respectively. The corresponding optimization problem is illustrated as

Object: fx

Subject: (1) fx0, (2) abfxdx=1, (3) fx1fx2=Mx1Mx2,x1,x2amin,bmax. Where subject (1) and (2) are the probability axiom that any PDF should obey and (3) is a constraint to balance the contribution for each interval, here Mx is a number function that x is comprised by M out of N intervals. To analytical solve this optimization problem, the following algorithm is developed.

Table 2. Algorithm for representation of pure interval uncertainty

Calculate the basic weight as Δ=1/i=1Nbi-ai
Count if xai,bi, for xamin,bmax, i=1:N
Compute the piecewise function Mx
Construct the converted PDF as fjx=ΔMjx,xamin,bmax, 1j2N-1

The converted PDF is actually a 2N-1 piecewise function and usually cannot be expressed in parametric form (e.g. Normal or Uniform distribution etc.).

5. Uncertainty quantification for a rotor/stator with mixed uncertainties

In the existing research works, most system parameters in the rotor are either set as fixed values [1-19] or assigned particular probability distributions [20-23]. However, in practical engineering, deterministic value or specific distribution is not always accessible due to various uncertainty factors. On the other hand, sparse data collected from previous/similar product and parameter intervals estimated based on expert judgment may be available. Taking full advantage of all these valuable information will lead an accurate analysis result

Considering these, three parameters with mixed uncertainty are taken into account. Specially, the rotation speed is set as a random variable which follows the Normal distribution as Ω~NμΩ,σΩ, while the available information for clearance R0 and friction coefficient f are given in a mixed sparse-data and interval form. We will concentrate on the effect of uncertainty for three parameters Ω, R0, f by setting other parameters fixed. From the uncertainty perspective, both aleatory uncertainty and epistemic uncertainty are taken into account, including the geometrical uncertainty in clearance, property uncertainty in friction coefficient and other uncertainties arising from manufacturing, mounting process etc.

For the sake of comparison, total 5 individual cases are demonstrated. They all have the same deterministic limit state function while the uncertain input variables are given in different styles. The developed likelihood-based approach is adopted to show its validity in uncertainty representation and quantification. The reliability of each case is calculated to quantify the effect to evaluation result of different format uncertain information. All available information is summarized in Table 3.

5.1. Case 1. Pure aleatory uncertainty

In this case, only pure aleatory uncertainty is discussed. The rotation speed Ω and clearance R0 are set as random variables, which are assumed to follow Normal distribution as Ω~NμΩ,σΩ and R0~NμR0,σR0. The marginal distribution for parameter Ω and R0 are plotted in Fig. 3 for reference and the joint PDF is constructed as Eq. (24):

p Ω , R 0 | Θ = p Ω | Θ p R 0 | Θ .

Table 3. Available information and data collection for case study

Available information
speed Ω
Aleatory (irreducible)
Complete Probability distribution
Ω ~ N ( μ Ω , σ Ω ) , μ= 0.8, σ= 0.04
friction coefficient f
Epistemic (physical)
Multiple input intervals
Source 1: [0.05, 0.1],
Source 2: [0.1, 0.15]
Source 3: [0.15, 0.2],
Source 4: [0.2, 0.3]
Clearance R0
Epistemic (geometry)
Probability distribution with uncertainty parameters and Multiple input intervals
R 0 ~ N ( μ , σ ) Source 1: [1.5, 1.8]
Source 2: [2, 2.2]
Source 3: [1.7, 2.1]
Sparse point data: 1.83, 1.91, 2.16

Fig. 3. Marginal distribution for parameter Ω and R0

 Marginal distribution for parameter Ω and R0

a) Rotation speed (Ω)

 Marginal distribution for parameter Ω and R0

b) Clearance (R0)

It is easily found that the 99.7 % credible intervals are within the area of [0.6, 1]×[1.5, 2.5], so this area is approximate to substitute for the infinite bound for the normal distribution.

As discussed in section 3, the rotor/stator system will lose its stability through saddle-node bifurcation boundary (SNu). Thus, the limit state function for this case is determined by the upper bound of saddle-node bifurcation as Eq. (25) and present in Fig. (4):

g ( R 0 , Ω , β , ζ , f ) = R 0 2 1 + f 2 - Ω 4 2 ζ Ω + f 2 + 1 + β - Ω 2 2
              - R 0 2 1 + β - Ω 2 + f 2 ζ Ω + f 2 .

The reliability is defined as the probability that this system works stable (in safe region) as Eq. (26). By generating 1×106 samples, it is calculated as 0.4195 through Monte-Carlo method. The safe/failure domain is intuitively depicted in Fig. 4:

R c a s e 1 = P s g ( R 0 , Ω , β , ζ , f ) > 0 .

In this scenario, only the stochastic information is considered, which is similar to previous research works considering uncertainties [20-22]. If other type uncertain information (e.g. interval) becomes available, the likelihood-based approach is adopted to perform an accurate reliability analysis.

5.2. Case 2. Pure epistemic uncertainty (single interval)

In this case, the effect to system response is investigated taking epistemic uncertainty into account. Both rotation speed Ω and clearance R0 are present in the interval form as Ω [0.6, 1], R0 [1.5, 2.5].

Fig. 4. Safe region and failure region for case 1

 Safe region and failure region for case 1

Fig. 5. Safe region and failure region for case 2

 Safe region and failure region for case 2

Since it is equal the chance of uncertain variable Ω and R0 that take any arbitrary value within the given intervals, the reliability can be expressed as the ratio of the area for safe region and that for whole interval region. The reliability is thus analytical derived as Eq. (27) and the safe region/failure region is also present in Fig. 5:

R c a s e 2 = A s A s + A f = Ω × R 0 g ( R 0 , Ω , β , ζ , f ) d Ω d R 0 l Ω × l R 0 = 0.4834 ,

where lΩ and lR0 correspond to the length of given interval for Ω and R0.

5.3. Case 3. Mixed uncertainties (single interval)

For this scenario, the rotation speed Ω is given as a random variable while the friction coefficient f is present as an interval variable. Specific, Ω~N(0.8,0.04) and f [0.1, 0.2].

Since no additional information is available for f, the uniform distribution is adopted to model the friction coefficient as f~U0.1,0.2. The joint PDF is constructed as Eq. (28) and present in Fig. 6:

p Ω , f | Θ = p Ω | Θ p f | Θ .

By setting other parameter fixed, the reliability is calculated as 0.4302 by drawing enough samples.

5.4. Case 4. Mixed uncertainties (multiple intervals)

In this case, the rotation speed Ω remains unchanged as a random variable Ω~N(0.8,0.04) while the friction f is given by multiple intervals: source 1: [0.05, 0.1], source 2: [0.1, 0.15], source 3: [0.15, 0.2], source 4: [0.2, 0.3].

By adopting the algorithm discussed in Section 4.4, the converted PDF is constructed as:

p = 10 / 3 , f 0.05,0.1 , 25 / 3 , f 0.1,0.15 , 5 , f 0.15,0.2 , 5 / 3 , f 0.2,0.3 ,

and shown in Fig. 7.

Fig. 6. Safe region and failure region for case 3

 Safe region and failure region for case 3

Fig. 7. Safe region and failure region for case 4

 Safe region and failure region for case 4

The reliability is finally obtained as 0.4165 by generating enough samples.

5.5. Case 5 Mixed uncertainties (multiple intervals with spares data)

In this case, the uncertainties in Ω, f and R0 are all taken into consideration. All available information listed in Table 3 is utilized to perform a credible reliability analysis.

The rotation speed Ω is given in stochastic form, thus no further transformation needs to take. To the uncertainty in f, the multiple input intervals from different sources can be dealt in similar way as case 4, so that samples can be easily generated from the PDF converted. In this case, we are focused on the uncertainty in R0.

Based on the accumulated knowledge, the uncertain parameter is assumed to follow one certain distribution (Normal) while the collected sparse data could be used to evaluate unknown statistic parameters (μ and σ). By adopting Bayesian rule, the likelihood of unknown parameters μ and σ is constructed as Eq. (29):

L Θ i = 1 m f X x i | Θ i = 1 n F X b i | Θ - F X b i | Θ ,

where fX and FX denote the PDF and CDF of normal distribution respectively. Θ is the parameter vector, for this case Θ=μ,σ, m=n= 3.

Since there is no additional available information, the non-informative prior is assigned as μ~Uniformaμ,bμ, σ~Uniformaσ,bσ, thus Θ) can be taken out of integration and the posterior joint PDF is reduced to:

p Θ | D = L Θ π Θ L Θ π Θ = i = 1 m f X x i | Θ i = 1 n F X b i | Θ - F X b i | Θ i = 1 m f X x i | Θ i = 1 n F X b i | Θ - F X b i | Θ d Θ .

There is usually no analytical solution for Eq. (30), so the posterior joint PDF is numerical evaluated. The marginal distributions of the parameter μ and σ are depicted in Fig. 8.

Fig. 8. PDF of mean and standard deviation

 PDF of mean and standard deviation

a) Mean of R0 (μ)

 PDF of mean and standard deviation

b) Standard deviation of R0 (σ)

Fig. 9. Unconditional PDF of R0

 Unconditional PDF of R0

Fig. 10. Limit state surface on the space of Ω, R0 and f

 Limit state surface on the space of Ω, R0 and f

The unconditional PDF is obtained by integrating over the space of parameters Θ as:

p R 0 | d a t a = p R 0 | Θ p Θ | d a t a d Θ .

Fig. (9) presents the unconditioned posterior PDF of variable R0.

The reliability is finally calculated as 0.6142 and the limit state surface is shown in the parameter space (Ω, f and R0) as Fig. 10.

It should be noted that even if the variable R0 is assumed to follow a particular type of a parametric distribution (Normal), the unconditional density after the integration in Eq. (31) (where the uncertainty in the estimates of the parameters Θ is accounted for) is non-parametric, i.e. the resultant probability distribution is not of the same type and cannot be classified as normal. See in Fig. 8.

6. Discussion

6.1. Benefits and achievements

In the context of uncertain dynamic analysis, the developed framework incorporates both deterministic data and uncertain information in a coherent way compared with different approaches for different kinds of data pursued in previous studies. This proposal is flexible and allows simultaneous processing of information in different formats, i.e. sparse point data, probability distribution and interval estimates, indicating a successful information aggregation and data fusion.

On the other hand, to achieve an accurate result, the system response function needs to be evaluated repeatedly to calculate the statistics of the output quantity. The uncertainties in both probability distribution and distribution parameters are included, which require a double loop sampling process. This procedure is cumbersome and may lead to extensive computation for large sample numbers. In the proposed approach, the computation cost is significantly reduced as the double loop sampling is collapsed into a single loop by discarding traditional Monte Carlo method. A unique PDF is constructed which is convenient for further dynamic analysis or risk assessment.

6.2. Validation and verification

It is found that the reliability analysis result differs significantly with different sources of uncertainty (case 1, 2 vs. case 3, 4) or representation of uncertainty in different formats (case 1 vs. case2, case 3 vs. case 4) taken into account, indicating a subtle influence of different uncertain variables to the outcome. Though these calculation results are derived based on numerical simulation, it still has profound significance in practical engineering. For example, the safe region obtained is smaller than expected when taking all possible uncertain factors into consideration. It is desired to increase the original clearance under deterministic situation for a safer rotor/stator design as rotating speed continues increasing. It should also be noticed that the calculation result (e.g. reliability) is only a reference value for an individual unit. Since the information used is given in statistical form (probability distribution), the result can only be experimentally validated under the circumstance with enough large samples.

7. Conclusions

In this paper, a modified Jeffcott rotor/stator system with mixed uncertainties is taken into study. The concept of likelihood is extended from classical point value to interval variables, which lay the foundation of proposed algorithm for PDF construction under mixed uncertainty. A likelihood-based approach addressing multiple inputs is developed, which enable a comprehensive uncertain dynamic analysis to the rotor/stator system incorporating all available information (sparse point data, probabilistic distribution and intervals). Several numerical cases are presented for demonstration and validation. Since the coexistence of both aleatory uncertainty and epistemic uncertainty is inevitable in practical engineering, the developed methodology has an extensive application prospect, such as industry manufacturing & mounting, risk assessment and product maintaining.


The authors gratefully acknowledge the financial support from the National Basic Research Program of China (973 Program, Project Number: 2013CB733000) and the National Natural Science Foundation of China (No. 51675026).


  1. Black H. F. Interaction of a whirling rotor with a vibrating stator across a clearance annulus. International Journal of Mechanical Engineering Science, Vol. 10, 1968, p. 1-12. [Publisher]
  2. Muszynska A. Synchronous self-excited rotor vibration caused by a full annular rub. Machinery Dynamics 8th Seminar, Halifax, Nova Scotia, Canada, Canadian Machinery Association, 1984, p. 1-21. [Search CrossRef]
  3. Muszynska A. Rotor-to-stationary element rub-related vibration phenomena in rotating machinery – literature survey. Sound and Vibration, Vol. 21, Issue 3, 1989, p. 3-11. [Search CrossRef]
  4. Bently D. E., Goldman P., Yu J. J. Full annular rub in mechanical seals. Part II: analytical study. Proceedings of ISROMAC-8, Honolulu, Hawaii, USA, 2000. [Search CrossRef]
  5. Jiang J., Ulbrich H. Stability analysis of sliding whirl in a nonlinear Jeffcott rotor with cross-coupling stiffness coefficients. Nonlinear Dynamics, Vol. 24, Issue 3, 2001, p. 269-283. [Publisher]
  6. Childs D. W. Rub-induced parametric excitation in rotors. Journal of Mechanical Design, Vol. 101, 1979, p. 640-644. [Publisher]
  7. Childs D. W. Fractional-frequency rotor motion due to non-symmetric clearance effects. Journal of Engineering for Power, Vol. 104, 1982, p. 533-541. [Publisher]
  8. Ehrich F. F. High order subharmonic response of high speed rotors in bearing clearance. ASME Journal of Vibration, Acoustics, Stress, and Reliability in Design, Vol. 110, 1988, p. 9-16. [Publisher]
  9. Day W. B. Asymptotic expansions in nonlinear rotor dynamics. Quarterly of Applied Mathematics, Vol. 44, Issue 4, 1987, p. 779-792. [Search CrossRef]
  10. Kim Y. B., Noah S. T. Quasi-periodic response and stability analysis for a nonlinear Jeffcott rotor. Journal of Sound and Vibration, Vol. 190, Issue 2, 1996, p. 239-253. [Publisher]
  11. Choi S. K., Noah S. T. Mode-locking and chaos in a Jeffcott rotor with bearing clearances. ASME Journal of Applied Mechanics, Vol. 61, 1994, p. 131-138. [Publisher]
  12. Goldman P., Muszynska A. Chaotic behavior of rotor/stator system with rubs. Journal of Engineering for Gas Turbins and Power, Vol. 116, 1994, p. 693-701. [Publisher]
  13. Chu F., Zhang Z. Bifurcation and chaos in a rub-impact Jeffcott rotor system. Journal of Sound and Vibration, Vol. 210, Issue 1, 1998, p. 1-18. [Publisher]
  14. Zhang W. Dynamic instability of multi-degree-of-freedom flexible rotor systems due to full annular rub. IMechE, 1988, p. 305-308. [Search CrossRef]
  15. Crandall S. From whirl to whip in rotor dynamics. Transactions of IFToMM 3rd International Conference on Rotor dynamics, Éd. du Centre NatiScient, 1990, p. 19-26. [Search CrossRef]
  16. Bartha A. R. Dry friction backward whirl of rotors: theory, experiments, results, and recommendations. 7th International Symposium on Magnetic Bearings, ETH Zurich, Schwitzland, 2000. [Search CrossRef]
  17. Bently D. E., Yu J. J., Goldman P. Full annular rub in mechanical seals. Part I: experimental results. Proceedings of ISROMAC-8, Honolulu, Hawaii, USA, 2000. [Search CrossRef]
  18. Jiang J., Ulbrich H. The physical reason and the analytical condition for the onset of dry whip in rotor-to-stator contact systems. Journal of Acoustics and Vibration – Transactions of the ASME, Vol. 127, 2005, p. 594-603. [Publisher]
  19. Jiang J. The analytical solution and existence condition of dry friction backward whirl in rotor-to-stator contact systems. Journal of Acoustics and Vibrations – Transactions of the ASME, Vol. 129, 2007, p. 260-264. [Publisher]
  20. Liao Haitao Global resonance optimization analysis of nonlinear mechanical systems: application to the uncertainty quantification problems in rotor dynamics. Communications in Nonlinear Science and Numerical Simulation, Vol. 19, Issue 9, 2014, p. 3323-3345. [Search CrossRef]
  21. Didier Jérome, Jacques Jean Sinou, Faverjon Béatrice Study of the non-linear dynamic response of a rotor system with faults and uncertainties. Journal of Sound and Vibration, Vol. 331, Issue 3, 2012, p. 671-703. [Search CrossRef]
  22. Jean-Jacques Sinou, Faverjon Béatrice The vibration signature of chordal cracks in a rotor system including uncertainties. Journal of Sound and Vibration, Vol. 331, Issue 1, 2012, p. 138-154. [Publisher]
  23. Motley Michael R., Young Yin L. Influence of uncertainties on the response and reliability of self-adaptive composite rotors. Composite Structures, Vol. 94, Issue 1, 2011, p. 114-120. [Publisher]
  24. Lingener A. Experimental investigation of reverse whirl of a flexible rotor. Transactions of IFToMM 3rd International Conference on Rotor dynamics, Lyon, France, 1990, p. 13-18. [Search CrossRef]
  25. Wu F., Flowers G. T. An experimental study of the influence of disk flexibility and rubbing on rotor dynamics. ASME, Vibration of Rotating Systems, Vol. 60, 1993, p. 19-26. [Search CrossRef]
  26. Helton J. C. Quantification of margins and uncertainties: conceptual and computational basis. Reliability Engineering and System Safety, Vol. 96, 2011, p. 976-1013. [Publisher]
  27. Helton J. C. Quantification of margins and uncertainties: alternative representations of epistemic uncertainty. Reliability Engineering and System Safety, Vol. 96, 2011, p. 1034-1052. [Publisher]
  28. Helton J. C., Johnson J. C., Sallaberry C. J. Quantification of margins and uncertainties: example analyses from reactor safety and radioactive waste disposal involving the separation of aleatory and epistemic uncertainty. Reliability Engineering and System Safety, Vol. 96, 2011, p. 1014-33. [Publisher]
  29. Zaman Kais, Rangavajhala Sirisha, McDonald Mark P., Mahadevan Sankaran A probabilistic approach for representation of interval uncertainty. Reliability Engineering and System Safety, Vol. 96, 2011, p. 117-130. [Publisher]
  30. Sankararaman Shankar, Mahadevan Sankaran Likelihood-based representation of epistemic uncertainty due to sparse point data and/or interval data. Reliability Engineering and System Safety, Vol. 96, 2011, p. 814-824. [Publisher]
  31. Wallstrom Timothy C. Quantification of margins and uncertainties: a probabilistic framework. Reliability Engineering and System Safety, Vol. 96, 2011, p. 1053-1062. [Publisher]
  32. Urbina Angel, Mahadevan Sankaran, Paez Thomas L. Quantification of margins and uncertainties of complex systems in the presence of aleatoric and epistemic uncertainty. Reliability Engineering and System Safety, Vol. 96, 2011, p. 1114-1125. [Publisher]
  33. Edwards A. W. F. Likelihood. 1st Ed. Cambridge University Press, Cambridge, 1972. [Search CrossRef]
  34. Pawitan Y. In All Likelihood: Statistical Modeling and Inference Using Likelihood. Oxford Science Publications, New York, 2001. [Search CrossRef]
  35. Barnard G. A., Jenkins G. M., Winsten C. B. Likelihood inference and time series. Journal of the Royal Statistical Society, Series A, Vol. 125, 1962, p. 321-372. [Publisher]

Cited By

Journal of Sound and Vibration
Chao Fu, Yuandong Xu, Yongfeng Yang, Kuan Lu, Fengshou Gu, Andrew Ball
Chao Fu, Dong Zhen, Yongfeng Yang, Fengshou Gu, Andrew Ball
Journal of Mechanical Science and Technology
Lechang Yang, Ketai He, Yanling Guo
Journal of Vibroengineering
Chao Fu, Xingmin Ren, Yongfeng Yang
Journal of Vibroengineering
Paweł Pietkiewicz, Sławomir Banaszek, Grzegorz Żywica