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

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:

The steady-state periodic solution is presented in the form as $\left\{\begin{array}{l}X=R\mathrm{c}\mathrm{o}\mathrm{s}\left(\mathrm{\Omega }\tau +\psi \right)\\ Y=R\mathrm{s}\mathrm{i}\mathrm{n}\left(\mathrm{\Omega }\tau +\psi \right)\end{array}\right.$. The amplitude and phase angle are given by:

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

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

By adopting the Routh-Hurwitz condition, it is found that the solution is always stable since $4{\zeta }^2\beta >0$. Thus, the two real roots ${\mathrm{\Omega }}_l$ and ${\mathrm{\Omega }}_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 $\mathbf{J}$ is time-dependent and is in following form:

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

where:

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

Define $\underset{\_}{\mathbf{J}}={\mathbf{T}}^{-1}\left(\mathbf{J}\cdot \mathbf{T}-d\mathbf{T}\right)$ and the $\mathbf{J}$ is solved as:

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

where $\lambda$ denotes the eigenvalues of $\mathbf{J}$ and coefficients are determined as:

As these eigenvalues are related to the amplitude $R\text{,}$ 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.

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 $p\left(x|\mathrm{\Theta }\right)$ where $\mathrm{\Theta }$ corresponds to the parameter set. The likelihood function $L\left(\mathrm{\Theta }|x\right)$ is defined as the probability of a series of given data set conditioned on distribution parameter $\mathrm{\Theta }$ [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 $x_i$ based on mean value theorem as Eq. (18) for practical use in previous work [30, 33, 34]:

The parameter $\mathrm{\Theta }$ 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 $\mathrm{\Theta }$ for interval [$a$, $b$] is:

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

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 $\mathrm{\Theta }$ 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 $\mathrm{\Theta }$ 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]:

Assigning an appropriate prior distribution, the parameter $\mathrm{\Theta }$ is estimated as Eq. (22) by adopting Bayesian theory:

where $\pi \left(\mathbf{\Theta }\right)$ denote the joint prior distribution of parameter vector $\mathbf{\Theta }$, while $p\left(\mathbf{\Theta }|D\right)$ is the estimated probability distribution of $\mathrm{\Theta }$. $l\left(\mathbf{\Theta }|D\right)$ is the combined likelihood function which has the same form as Eq. (21).

It should be noted that the parameter $\mathrm{\Theta }$ 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 $\theta$. 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 $\theta$ are generated in the outer loop, each of which corresponds to a distribution of $p\left(\theta \right)$, and $N$ samples of $X$ are drawn for each sample $\theta$ in the inner loop, the computation cost is total $M$×$N$. To get a satisfactory result, 10⁴ samples is usually required for $M$ and $N$, thus the number of total samples is incredible large.

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 $p\left(X|\theta \right)$ with a unique distribution $p\left(X\right)$:

The unconditional PDF $p\left(X\right)$ can be calculated by integrating the conditional PDF $p\left(x|\mathbf{\Theta }\right)$ over the entire parameter space of $\mathrm{\Theta }$.

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. $X\in \left[a,b\right]$. From the probabilistic perspective, it can be interpreted as $X~U\left(a,b\right)$ 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: $f\left(x\right)$

Subject: (1) $f\left(x\right)\ge 0$, (2) $\underset{a}{\overset{b}{\int }}fr\left(x\right)dx=1$, (3) $f\left(x\right)\equiv p$.

The solved PDF $f\left(x\right)$ is the corresponding probabilistic description for the given interval [$a$, $b$].

For multiple $N$ input intervals $x\in \left[a_i,b_i\right]$, $1\le i\le N$, there are total $2N-1$ subintervals $\left[d_j,d_{j+1}\right]$, $a_{min}\le d_j\le d_{j+1}\le b_{max}$, where $a_{min}$ is the lower bound and $b_{max}$ is the upper bound respectively. The corresponding optimization problem is illustrated as

Object: $f\left(x\right)$

Subject: (1) $f\left(x\right)\ge 0$, (2) $\underset{a}{\overset{b}{\int }}f\left(x\right)dx=1$, (3) $\frac{f\left(x_1\right)}{f\left(x_2\right)}=\frac{M\left(x_1\right)}{M\left(x_2\right)},\forall x_1,x_2\in \left[a_{min},b_{max}\right]$. 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 $M\left(x\right)$ 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.

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.1. Case 1. Pure aleatory uncertainty

In this case, only pure aleatory uncertainty is discussed. The rotation speed $\mathrm{\Omega }$ and clearance $R_0$ are set as random variables, which are assumed to follow Normal distribution as $\Omega ~N\left({\mu }_{\mathrm{\Omega }},{\sigma }_{\mathrm{\Omega }}\right)$ and $R_0~N\left({\mu }_{R_0},{\sigma }_{R_0}\right)$. The marginal distribution for parameter $\mathrm{\Omega }$ and $R_0$ are plotted in Fig. 3 for reference and the joint PDF is constructed as Eq. (24):

Fig. 3Marginal distribution for parameter Ω and R0

a) Rotation speed ($\mathrm{\Omega }$)

b) Clearance ($R_0$)

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 ($S{N}_u$). 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):

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

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 $\mathrm{\Omega }$ and clearance $R_0$ are present in the interval form as $\mathrm{\Omega }\in$ [0.6, 1], $R_0\in$[1.5, 2.5].

Fig. 4Safe region and failure region for case 1

Fig. 5Safe region and failure region for case 2

Since it is equal the chance of uncertain variable $\mathrm{\Omega }$ and $R_0$ 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:

where $l_{\mathrm{\Omega }}$ and $l_{R0}$ correspond to the length of given interval for $\mathrm{\Omega }$ and $R_0$.

5.3. Case 3. Mixed uncertainties (single interval)

For this scenario, the rotation speed $\mathrm{\Omega }$ is given as a random variable while the friction coefficient $f$ is present as an interval variable. Specific, $\mathrm{\Omega }~N\left(\mathrm{0.8,0.04}\right)$ and $f\in$ [0.1, 0.2].

Since no additional information is available for $f$, the uniform distribution is adopted to model the friction coefficient as $f~U\left(\mathrm{0.1,0.2}\right)$. The joint PDF is constructed as Eq. (28) and present in Fig. 6:

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 $\mathrm{\Omega }$ remains unchanged as a random variable $\mathrm{\Omega }~N\left(\mathrm{0.8,0.04}\right)$ 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:

and shown in Fig. 7.

Fig. 6Safe region and failure region for case 3

Fig. 7Safe 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 $\mathrm{\Omega }$, $f$ and $R_0$ are all taken into consideration. All available information listed in Table 3 is utilized to perform a credible reliability analysis.

The rotation speed $\mathrm{\Omega }$ 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 $R_0$.

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 ($\mu$ and $\sigma$). By adopting Bayesian rule, the likelihood of unknown parameters $\mu$ and $\sigma$ is constructed as Eq. (29):

where $f_X\left(\cdot \right)$ and $F_X\left(\cdot \right)$ denote the PDF and CDF of normal distribution respectively. $\mathbf{\Theta }$ is the parameter vector, for this case $\mathbf{\Theta }=\left(\mu ,\sigma \right)$, $m=n=$ 3.

Since there is no additional available information, the non-informative prior is assigned as $\mu ~Uniform\left(a_{\mu },b_{\mu }\right)$, $\sigma ~Uniform\left(a_{\sigma },b_{\sigma }\right)$, thus $\mathbf{\Theta }$) can be taken out of integration and the posterior joint PDF is reduced to:

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

Fig. 8PDF of mean and standard deviation

a) Mean of $R_0$ ($\mu$)

b) Standard deviation of $R_0$ ($\sigma$)

Fig. 9Unconditional PDF of R0

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

The unconditional PDF is obtained by integrating over the space of parameters $\mathrm{\Theta }$ as:

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

The reliability is finally calculated as 0.6142 and the limit state surface is shown in the parameter space ($\mathrm{\Omega }$, $f$ and $R_0$) as Fig. 10.

It should be noted that even if the variable $R_0$ 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 $\mathrm{\Theta }$ 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.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.