Steady-State Response Analysis of an Uncertain Rotor Based on Chebyshev Orthogonal Polynomials

Abstract

The performance of a rotor system is influenced by various design parameters that are neither precise nor constant. Uncertainties in rotor operation arise from factors such as assembly errors, material defects, and wear. To obtain more reliable analytical results, it is essential to consider these uncertainties when evaluating rotor performance. In this paper, the Chebyshev interval method is employed to quantify the uncertainty in the steady-state response of the rotor system. To address the challenges of high-dimensional integration, an innovative sparse-grid integration method is introduced and demonstrated using a rotor tester. The effects of support stiffness, mass imbalance, and uncertainties in the installation phase angle on the steady-state response of the rotor system are analyzed individually, along with a comprehensive assessment of their combined effects. When compared to the Monte Carlo simulation (MCS) method and the full tensor product grid (FTG) method, the proposed method requires only 68% of the computational cost associated with MCS, while maintaining calculation accuracy. Additionally, sparse-grid integration reduces the computational cost by approximately 95.87% compared to the FTG method.

1. Introduction

The rotor system is a vital component of rotating machinery such as aero engines and gas turbines. In the field of rotor dynamics, numerical simulators are commonly used to predict the dynamic behavior of the rotor system. However, uncertainties inherent in actual rotor systems stem from operational conditions and geometric parameters that may vary, including loose connections [1] or cracks [2]. These variations can lead to deviations from the expected steady-state response. Consequently, a steady-state response of the rotor system which initially meets design specifications may exceed acceptable limits in practice. A major challenge for engineers and researchers is to analyze and quantify the impacts of these uncertainties on the steady-state responses of complex rotor systems [3,4,5].

There are four commonly used methods for uncertainty analysis: the stochastic finite element method, probabilistic models, fuzzy methods, and interval methods. The stochastic finite element method assesses uncertainty by establishing stochastic differential equations for the system’s response based on stochastic matrix theory or Taylor’s formula [6]. Kureishi [7] applied this method to analyze a multi-degree-of-freedom rotor system in both time and frequency domains, investigating dynamics such as axial trajectory and random amplitude and frequency responses. Liu [8] developed a stochastic finite element model of a rotor-bearing system considering radial clearance randomness and employed the orthogonal polynomial method to examine the effects of factors like mean and coefficient of variation of radial clearances, eccentricity, and rotational speeds on the system’s dynamic performance. However, the stochastic finite element method’s disadvantage is found in the difficulty and time consumed in deriving embedded equations when analyzing complex rotor systems. Probabilistic methods necessitate complete probability distributions of uncertain parameters, which are often challenging to obtain [9]. Fuzzy methods employ fuzzy variables to represent uncertain parameters, offering a robust framework for capturing the inherent ambiguity in a system’s characteristics. The Fuzzy Stochastic Finite Element Method (FSFEM) is particularly advantageous, as it integrates fuzzy logic with stochastic modeling, allowing for a more comprehensive analysis of uncertainties in dynamic systems. Lara [10,11] utilized the FSFEM to investigate flexible rotor dynamics under uncertain parameters, examining the impacts of variables such as pad radius, oil viscosity, and radial clearance on tilting-pad journal bearings. However, the derivation of membership functions in fuzzy methods can be complex, presenting challenges when applied to the intricate task of rotor uncertainty analysis. Additionally, the effectiveness of the FSFEM may be limited by the presence of multi-source uncertainties, which complicates the modeling process and can hinder accurate assessments of system behavior.

Interval analysis is widely applied in uncertainty analysis. Although variations in system parameters are typically random, their range is usually finite, described as “unknown but bounded”, indicating that the parameters fluctuate within definable intervals [12]. In the interval analysis of structural dynamics, uncertain parameters are represented as interval vectors, and differential equations governing rotor motion are formulated by incorporating these interval vectors. The solutions to these equations are derived using interval mathematics and regression theory, which facilitates the determination of the range of the system’s steady-state response [13]. Ma [14,15] proposed a Bound Correction Interval Analysis Method (BCIM) based on Chebyshev expansion and modal superposition for the analysis of uncertain rotors. Mao [16] combined interval regression analysis with regularization to introduce a computational inverse method for identifying bearing loads. Utilizing interval analysis methods, Fu [17,18,19] examined the effects of uncertainties in several design parameters on rotors, including stiffness, damping, Young’s modulus, cracking, and coupling misalignment. The interval analysis approach allows for the omission of specific probability distribution models for the uncertain parameters, requiring only the definition of their upper and lower fluctuation bounds. However, in scenarios characterized by multi-source uncertainty, this method necessitates the computation of combinations of specified integration points, such as Gaussian integration points, for multi-dimensional interval variables, thus significantly increasing the computational workload. To enhance computational efficiency, dimensionality reduction can be achieved through sparse generalized polynomial chaos expansion (gPCE) or polynomial dimensional decomposition [20,21,22]. Lu [23] employed a polynomial dimensional decomposition method to analyze the uncertainty in the response of a nonlinear rotor system, achieving effective dimensionality reduction. Nonetheless, when considering the nonlinear factors of the rotor system (e.g., nonlinear oil film damping), balancing computational efficiency and accuracy remains an area requiring further investigation.

The remainder of this paper is organized as follows: Section 2 presents the proposed methodology; this is followed by Section 3, which introduces an example of the research subject, namely a dynamically similar rotor tester for a turboshaft engine. Section 4 then discusses the results of the numerical analysis. Finally, Section 5 concludes the paper.

2. The Proposed Method

This section describes the method proposed in this paper, which encompasses both Chebyshev orthogonal polynomial approximation and the fundamental techniques of sparse-grid numerical integration. In the context of high-dimensional random variables, the use of sparse-grid numerical integration can effectively reduce computation time while maintaining calculation accuracy, thereby alleviating the “curse of dimensionality”. To facilitate a clearer understanding, a flowchart is provided below for conducting uncertainty analysis of the steady-state response of rotor systems based on Chebyshev orthogonal polynomials.

2.1. Chebyshev Orthogonal Polynomial Approximation

Suppose $\xi =[{\xi }_1,{\xi }_2,\cdots ,{\xi }_n]$, ξ_j (j = 1, $\cdots$, n) is a standard random variable with a uniform distribution $\mathcal{U}(-1, 1)$. Using the tensor product approach, an n-dimensional, k-order Chebyshev polynomial $T_{i_1,\cdots ,i_n}\left(\xi \right)$ can be defined as follows:

$$T_{i_1,\cdots ,i_n}\left(\xi \right)={\displaystyle \prod _{i_1+\cdots +i_n=k}{T}_{i_1}\left({\xi }_1\right)\cdots T_{i_n}\left({\xi }_n\right)}$$

(1)

where i₁, i₂, …, i_n are non-negative integers. $T_{i_j}$ denotes a Chebyshev polynomial with order i_j. For example, if n = 2, then $\xi =\{{\xi }_1,{\xi }_2\}$, and we have

$$\begin{array}{ll}k=0:& T_{0,0}(\xi )=T_0({\xi }_1)T_0({\xi }_2)=1\\ k=1:& T_{1,0}(\xi )=T_1({\xi }_1)T_0({\xi }_2)={\xi }_1, L_{0,1}(\xi )=T_0({\xi }_1)T_1({\xi }_2)={\xi }_2,\\ k=2:& T_{2,0}(\xi )=T_2({\xi }_1)T_0({\xi }_2)=2{\xi }_1^2-1, T_{1,1}(\xi )=T_1({\xi }_1)T_1({\xi }_2)={\xi }_1{\xi }_2,\\ & T_{0,2}(\xi )=T_0({\xi }_1)T_2({\xi }_2)=2{\xi }_2^2-1\\ \cdots \end{array}$$

(2)

A Chebyshev orthogonal polynomial of order p is used to approximate the multi-dimensional continuous function $f\left({\xi }_1,\cdots ,{\xi }_n\right):{\Re }^n\to \Re$, which can be expressed as follows:

$$\begin{array}{ll}f\left(\xi \right)\approx & {\displaystyle \sum _{i_1+\cdots +i_n\le p}{c}_{i_1,\cdots ,i_n}{T}_{i_1,\cdots ,i_n}\left(\xi \right)}={\displaystyle \sum _{i_1+\cdots +i_n=0}{c}_{i_1,\cdots ,i_n}{T}_{i_1,\cdots ,i_n}\left(\xi \right)}+\\ & {\displaystyle \sum _{i_1+\cdots +i_n=1}{c}_{i_1,\cdots ,i_n}{T}_{i_1,\cdots ,i_n}\left(\xi \right)}+\cdots +{\displaystyle \sum _{i_1+\cdots +i_n=p}{c}_{i_1,\cdots ,i_n}{T}_{i_1,\cdots ,i_n}\left(\xi \right)}\end{array}$$

(3)

where N_c is the total number of undetermined coefficients, which is given by the following equation:

$$N_c=1+{\displaystyle \frac{n!}{(n-1)!}}+{\displaystyle \frac{\left(n+1\right)!}{\left(n-1\right)!2!}}+{\displaystyle \frac{\left(n+2\right)!}{\left(n-1\right)!3!}}+\cdots +{\displaystyle \frac{\left(n-1+p\right)!}{\left(n-1\right)!p!}}={\displaystyle \frac{\left(n+p\right)!}{n!p!}}$$

(4)

This weight function serves to guarantee that Chebyshev polynomials of varying orders remain orthogonal over the specified interval. For the standard Chebyshev polynomials, the weight function is expressed as follows:

$$\mathcal{W}\left(\xi \right)={\displaystyle \frac{1}{\sqrt{1-{\xi }_1^2}\cdots \sqrt{1-{\xi }_n^2}}}$$

(5)

Multiplying both sides of Equation (3) by $T_{i_1,\cdots ,i_n}\left(\xi \right)$ and integrating the weight function $\mathcal{W}\left(\xi \right)$, the following results can be obtained:

$${\displaystyle {\int }_{-1}^1\cdots {\displaystyle {\int }_{-1}^1{T}_I\left(\xi \right)T_J\left(\xi \right)\mathcal{W}\left(\xi \right)\mathrm{d}{\xi }_1\cdots \mathrm{d}{\xi }_n}}=\left\{\begin{array}{cc}{\left({\displaystyle \frac{1}{2}}\right)}^l{\pi }^n& I=J\\ 0& I\ne J\end{array}\right.$$

(6)

where l is the number of non-zeroes in $i_1,\cdots ,i_n$. The undetermined coefficients in Equation (3) can be formulated as

$$c_{i_1,\cdots ,i_n}={\displaystyle \frac{2^l}{{\pi }^n}}{\displaystyle {\int }_{-1}^1\cdots {\displaystyle {\int }_{-1}^{1}f\left(\xi \right)T_{i_1,\cdots ,i_n}\left(\xi \right)\mathcal{W}\left(\xi \right)\mathrm{d}{\xi }_1\cdots \mathrm{d}{\xi }_n}}$$

(7)

When the dimension n is small (n ≤ 5), the coefficient in Equation (7) can be calculated by full factor numerical integration, and the calculation formula is as follows:

$$c_{i_1,\cdots ,i_n}\approx {\displaystyle \frac{2^l}{{\pi }^n}}{\displaystyle \sum _{j_1=1}^q\cdots {\displaystyle \sum _{j_n=1}^{q}f\left({\xi }_1^{\left(j_1\right)},\cdots ,{\xi }_n^{\left(j_n\right)}\right)}}\left(w_1^{\left(j_1\right)}\cdots w_n^{\left(j_n\right)}\right)$$

(8)

where ${\xi }_1^{\left(j_1\right)},\cdots ,{\xi }_n^{\left(j_n\right)}$ is integration node, and q is the number of integration nodes in each dimension.

2.2. Sparse Grid Numerical Integration

In Equation (8), if the number of nodes taken for each dimensional random variable is the same and equal to q, then the total number of function evaluations required is qⁿ. This approach for calculating $c_{i_1,\cdots ,i_n}$ is referred to as full-tensor-product grid (FTG). Regarding the selection of the number of nodes q (usually q > 3) in each dimension, a larger q results in higher accuracy; however, it also leads to an increase in computational workload. Consequently, when the dimension n exceeds 5, the strategy encounters a significant issue known as the “curse of dimensionality”. To alleviate this issue, an efficient approach to computing high-dimensional integrals is to use numerical integration with sparse-grid (SG) technology. The sparse grid numerical integration, also known as the sparse tensor product grid, was first proposed by Smolyak [24] and has been widely used by researchers to address high-dimensional integration problems, thereby improving the efficiency of integration. The SG with accuracy level L and n-dimensional integration is denoted as $\mathcal{G}(L,n)$, and its construction formula is as follows:

$$\mathcal{G}(L,n)={\displaystyle \underset{\max (n,L+1)\le {‖k‖}_1\le L+n}{\cup }\underset{i=1}{\overset{n}{\otimes }}{U}_1^{k_i}}$$

(9)

where ${‖k‖}_1={\displaystyle {\sum }_{i=1}^n{k}_i}$ (k_i ≥ 1) is the 1-norm of k, and k_i is a positive integer that determines the number of integration nodes assigned to each dimensional random variable. The accuracy level L is defined as q = L + 1. $U_1^{k_i}$ represents a univariate quadrature rule with k_i integral nodes. The SG decomposes the direct tensor product into different sub-tensor products, using only the sub-tensor products corresponding to multi-indices that satisfy max(n, L + 1) ≤ ${‖k‖}_1$ ≤ L + n. For example, if n = 2 and L = 2, the required multi-indices are k = {(1,2), (1,3), (2,1), (2,2), (3,1)}, and the construction of the corresponding sparse grid $\mathcal{G}(2,2)$ is given in Figure 1.

In Figure 1, the dots represent integration nodes; each square box represents a [−1,1]² space. The Gauss–Chebyshev integration points are used for each one-dimensional integration node, and the sub-grid corresponding to the multi-exponential k that satisfies the requirement is between the two dot–dash lines. Note that the tensor product results for sub-grids $U_1^1\otimes U_2^3$ and $U_1^3\otimes U_2^1$ in Figure 1a have one identical integration node, and the corresponding number of FTG nodes is 9. However, the number of integral nodes of $\mathcal{G}(2,2)$ in Figure 1b is 13. In lower-dimensional cases (n = 2), achieving the same integration accuracy (L = 2) requires more integration nodes for the SG method, compared to the FTG method (13 > 9). This discrepancy arises because SG techniques are specifically designed for high-dimensional problems (typically n > 5).

Table 1. Comparison of the number of integration nodes for SG and FTG.

n	L	SG	FTG	L	SG	FTG
5	1	11	32	2	66	243
7	1	15	128	2	113	2187
10	1	21	1024	2	231	59,049
13	1	27	8192	2	243	1,594,323

presents the number of integration nodes required for both SG and FTG methods, corresponding to different dimensions and accuracy requirements.

Table 1 demonstrates that as dimensionality increases, the number of integration nodes required for the FTG method rises rapidly, while the growth rate of required nodes in the SG method is significantly slower. This results in a substantial reduction in the number of integration nodes compared to the full-grid method, effectively alleviating the “curse of dimensionality”.

SG numerical integration is constructed from tensor products of hierarchical difference sets based on a one-dimensional product rule. The formula for integration based on the constructed sparse grid $\mathcal{G}(L,n)$ is as follows:

$$\mathcal{G}(L,n)\left[f\right]={\displaystyle \sum _{\max (n,L+1)\le {‖\mathit{k}‖}_1\le L+n}b({‖\mathit{k}‖}_1)\cdot (U_1^{k_1}\otimes \cdots \otimes U_1^{k_n})}\left[f\right]$$

(10)

where f is an n-dimensional (n > 1) continuous function, $U_1^{k_i}$($U_1^0=0$) is one-dimensional product rule, and coefficient $b({‖\mathit{k}‖}_1)$ is given by

$$b({‖\mathit{k}‖}_1){(-1)}^{L+n-{‖k‖}_1}\left(\begin{array}{c}n-1\\ L+n-{‖\mathit{k}‖}_1\end{array}\right)$$

(11)

Assuming that the number of n-dimensional SG integration nodes is N, denoted as $\mathcal{G}(L,n)$ = {l₁, …, l_i, …, l_N}, and the weight corresponding to the ith n-dimensional integration node l_i is represented as w_i, Equation (10) can be simplified to the following form:

$$\mathcal{G}(L,n)\left[f\right]={\displaystyle \sum _{i=1}^{N}f\left({\mathit{l}}_i\right)\cdot w_i}$$

(12)

The weight w_i is calculated by the following equation:

$$w_i={(-1)}^{L+n-{‖k_i‖}_1}\left(\begin{array}{c}n-1\\ L+n-{‖k_i‖}_1\end{array}\right)(w_1^{\left(i\right)}\cdots w_n^{\left(i\right)})$$

(13)

where $w_1^{\left(i\right)}\cdots w_n^{\left(i\right)}$ represents the product of the corresponding weights of each one-dimensional integration node in the ith dimensional sparse-grid integration node.

2.3. Uncertainty Analysis Flow Based on Chebyshev Orthogonal Polynomials

The process of interval uncertainty analysis using Chebyshev orthogonal polynomials consists of four primary steps. First, the random variables are specified, and one-dimensional interpolation points are generated using Gauss integration nodes. These nodes serve as the foundation for the subsequent steps in the analysis. Second, n-dimensional interpolation points are created from each one-dimensional point using either the tensor product or sparse-grid integration method. The choice between these two methods depends on the complexity of the problem and the desired balance between accuracy and computational efficiency. Once the n-dimensional interpolation points are established, the system response is calculated at each point, providing a comprehensive set of data for the uncertainty analysis.

In the third step, the coefficients of the Chebyshev orthogonal polynomials are computed to construct an orthogonal polynomial approximation model. This model serves as a surrogate for the original system, allowing for a more efficient uncertainty analysis. Finally, the steady-state response bounds of the rotor systems are calculated using the constructed approximation model. Once the Chebyshev orthogonal polynomial approximation is established, it can be used as a proxy model for the system, enabling efficient and accurate estimation of the steady-state response bounds through random sampling of the standard random variables using the MCS approach. Compared to the original finite element model, the use of Chebyshev orthogonal polynomials significantly enhances the efficiency of uncertainty analysis by reducing the computational burden while maintaining a high level of accuracy. Figure 2 provides a visual representation of this procedure, illustrating the flow of data and the key steps involved in the interval uncertainty analysis using Chebyshev orthogonal polynomials, thereby facilitating a clearer understanding of the methodology employed.

The accuracy and predictive quality of the Chebyshev orthogonal polynomial surrogate model can be assessed using the relative generalization error, denoted as ε_r. The formula for calculating ε_r is as follows:

$${\epsilon }_r={\displaystyle \frac{\mathrm{E}\left[{\left(\mathcal{M}({\mathcal{X}}_{\mathrm{val}})-{\mathcal{M}}_{Cheby}({\mathcal{X}}_{\mathrm{val}})\right)}^2\right]}{\mathrm{Var}\left[\mathcal{M}({\mathcal{X}}_{\mathrm{val}})\right]}}$$

(14)

where E[∙] denotes the expectation operator, while Var[∙] represents the variance operator. $\mathcal{M}$ represents the finite element model of the rotor system (also referred to as the original model or source model), ${\mathcal{M}}_{Cheby}$ denotes the Chebyshev orthogonal polynomial, and ${\mathcal{X}}_{\mathrm{val}}$ represents the random validation set.

3. Rotor System

Figure 3 illustrates the dynamics similarity tester for a turboshaft engine rotor with the discs and shafts constructed from the same material (40CrMnSiA), which has an elastic modulus of 2.11 × 10¹¹ N/m². The shaft has a length of 2.4 m and a diameter of 0.17 m. The masses of Disc 1, Disc 2, and Disc 3 are 3.3 kg, 5.83 kg, and 3.09 kg, respectively, with corresponding moments of inertia of 0.0079 kg·m², 0.0247 kg·m², and 0.0088 kg·m². The damping coefficient is 200 Ns/m. Other relevant model parameters include support stiffness, unbalanced mass on each disc, and installation phase angles. The installation phase angle of the unbalanced mass denotes the angular positions of each disc at the time of installation. The arrangement of these phase angles directly impacts the relative positioning of the unbalanced masses, thereby influencing the system’s vibration patterns and stability. Specific values and interpretations of these parameters in the example rotor system are presented in

Table 2. Model parameters.

Parameters	Dimension	Description
k1	2.46 × 105 N/m	Stiffness of Support 1
k2	3.64 × 105 N/m	Stiffness of Support 2
u1	30 g·cm	Mass unbalance on Disc 1
u2	30 g·cm	Mass unbalance on Disc 2
u3	30 g·cm	Mass unbalance on Disc 3
ϕ1	0°	Phase angle of unbalanced mass on Disc 2 relative to Disc 1
ϕ2	0°	Phase angle of unbalanced mass on Disc 3 relative to Disc 1

.

The steady-state response of the rotor system is computed using the finite element method (FEM). For clarity, the detailed modeling process is presented in Appendix A, ensuring that the organization of this paper remains coherent. Using the established FEM, the first three critical speeds are identified as 1434 rpm, 1964 rpm, and 6923 rpm. The corresponding vibration mode shapes are depicted in Figure 4, in which the first and second modes represent pitching, while the third mode represents bending. Taking the steady-state response associated with Disc 1 as an example, the frequency response function for Disc 1 is obtained, as illustrated in Figure 5.

4. Numerical Results

In this section, the effects of various uncertainty parameters, including support stiffness, unbalance, and phase, as well as their combinations, are investigated with respect to their relationship to the steady-state response of the rotor system, using the method proposed in this paper. The range of response variation is obtained and compared with the results calculated through the MCS method.

4.1. Effect of Interval Support Stiffness

Assuming the degree of uncertainty is ${\beta }_{K_b}$, the support stiffness can be expressed in interval form as follows:

$$\left[{\underset{\_}{K}}_b,{\overline{K}}_b\right]=\left[K_b-{\beta }_{K_b}{K}_b,K_b+{\beta }_{K_b}{K}_b\right]$$

(15)

where $K_b$ is the nominal value of stiffness in Support 1 or Support 2.

Generally, support stiffness comprises three components: the support structure stiffness, the squeeze film damper stiffness, and the bearing stiffness; all of these vary with temperatures, assembly states, and loads [25]. Furthermore, squirrel-cage elastic supports are widely used in jet engines, and unexpected breakage of the cage bars further increases the uncertainty of the support stiffness. Consequently, defining the precise value of support stiffness is challenging. Here, uncertain degrees of 5%, 10%, and 15% for support stiffness are considered, and uncertainty analysis is performed for these three scenarios to determine the range of steady-state responses of the rotor system. The convergence results of the relative generalization error ε_r are shown in Figure 6.

As shown in Figure 6, the error ε_r converges when the Chebyshev orthogonal polynomial is of order 3 (p = 3). Therefore, a third-order polynomial with four integration points in each dimension is employed. Utilizing full factorial numerical integration requires 16 invocations of the system evaluation model. The results of the uncertainty analysis for the steady-state response of Disc 1 are presented in Figure 7.

As shown in Figure 7, near the critical speed, the range of vibration amplitude broadens as the degrees of uncertainty increase. This phenomenon is more pronounced at the first and second critical speeds compared to the third, indicating that the dynamic response at these speeds is more sensitive to changes in support stiffness. In the range from 1000 to 4000 rpm, the resonance frequency, initially characterized by a single-peak mode, transitions into a broad resonance band. At the upper bound, the distinction between the first and second order becomes very blurred. The steady-state response exhibits less variability when the rotational speed exceeds 4000 rpm.

4.2. Effect of Interval Mass Unbalance

Unbalance in a rotor is unavoidable, due to manufacturing and assembly errors. The mass unbalance is considered as an interval variable, one expressed in the following form:

$$\left[\underset{\_}{e},\overline{e}\right]=\left[e-{\beta }_{e}e,e+{\beta }_{e}e\right]$$

(16)

where e is the nominal value (30 g·cm) of mass unbalance in Discs 1, 2, and 3, and ${\beta }_e$ is the degree of uncertainty. The steady-state response of the rotor system is calculated for uncertainty degrees of 5%, 10%, and 15%, respectively. Similar to the convergence analysis described in Section 4.1, a third-order (p = 3) Chebyshev orthogonal polynomial with four Gaussian integration nodes in each dimension has been constructed. The variation in the steady-state response of Disc 1 is shown in Figure 8.

Figure 8 illustrates that the range of steady-state amplitude widens as the rotational speed increases, indicating that higher rotational speeds amplify the effect of mass unbalance uncertainty on the vibration amplitude of the rotor system. Additionally, the upper and lower bounds are symmetric about the deterministic curve, while the critical speed remains unchanged.

4.3. Effect of Interval Unbalance Phase

Taking the unbalance phase of Disc 1 as a reference and assuming that the unbalance phase of Disc 1 is 0°, the unbalance phase shifts of Discs 2 and 3 fluctuate within the range of [0, 360°]. The convergence results of the relative error at different speeds are presented in Figure 9.

As shown in Figure 9, to consider the convergence of the relative error, a Chebyshev orthogonal polynomial of order (p = 8) is constructed. This forms the basis for conducting an uncertainty propagation analysis of the phase angle of the mass imbalance across the entire speed range. Figure 10 illustrates the variation range of the steady-state response at Disc 1.

As shown in Figure 10, uncertainty in the unbalance phase significantly impacts the steady-state response of the rotor, particularly at the third-order critical speed, at which the maximum amplitude increases by nearly threefold compared to the deterministic case. Changes in the unbalance phase alter the direction of the corresponding centrifugal force, consequently modifying the total centrifugal force on the rotor and the force couples. To ensure that the rotor operates reliably and safely, it is essential to account for the uncertainty of the unbalance phase during the design and assembly process, thereby preventing excessive unbalance-phase differences between the discs due to oversight.

4.4. Effect of Multi Interval Parameters

In practice, the uncertain parameters affecting rotor performance are typically diverse. This subsection investigates the simultaneous effects of seven uncertainty parameters—namely, interval stiffness in supports 1 and 2, interval mass unbalance in Discs 1, 2, and 3, and interval unbalance phase in Discs 2 and 3—on the rotor’s steady-state response. The degrees of uncertainty for stiffness and mass unbalance are assumed to be 5%, 10%, and 15%, respectively. Given the significant impact of phase fluctuations on vibration, only the scenario where the phase angle varies within ±5° is considered, with the phase of Disc 1 fixed as a reference. Three cases are established based on the degrees of uncertainty for all parameters, and uncertainty analysis is conducted for these scenarios. For clarity, the degrees of uncertainty for all parameters are summarized in

Table 3. The degrees of uncertainty for all parameters.

Parameters	Degrees of Uncertainty
	Case 1	Case 2	Case 3
k1	5%	10%	15%
k2	5%	10%	15%
u1	5%	10%	15%
u2	5%	10%	15%
u3	5%	10%	15%
ϕ1	±5°	±5°	±5°
ϕ2	±5°	±5°	±5°

.

When the Chebyshev orthogonal polynomial is of order 3 (p = 3), the relative error ε_r achieves convergence. Consequently, a third-order Chebyshev orthogonal polynomial is utilized to calculate the range of steady-state response variations of the system under multiple sources of uncertainty and varying parameter intervals. The range of steady-state response for Disc 1, obtained using the proposed method, is illustrated in Figure 11.

As shown in Figure 11, having multiple uncertainty parameters significantly affects the steady-state response of the rotor. This pattern of influence encompasses the characteristics exhibited by the individual parameters described earlier, including broad resonance bands, peak shifts, and fluctuations in amplitude magnitude.

4.5. Comparison and Validation

To validate the correctness of the proposed method, a comparison is made using the MCS method. Since the comparison processes for different uncertainties are consistent, due to space constraints, only one scenario is analyzed here, specifically, the one in which the uncertainties in support stiffness and mass imbalance are set at 10%, while the phase angles in Discs 2 and 3 fluctuate within a range of ±5°. The comparison’s results are illustrated in Figure 12 and Figure 13.

In Figure 12, the random sample size used for the MCS is 1000, and L in Figure 13 is the accuracy level of the sparse-grid integration used for establishing the Chebyshev orthogonal polynomial. The results indicate that the envelope intervals obtained by the two methods coincide, except at the anti-resonance point near 6500 rpm. Additionally, the statistical means and standard deviations of the steady-state responses obtained by both methods are nearly identical, which verifies the accuracy of the method proposed in this paper under conditions of multiple uncertainties.

Table 4. Comparison of computational efficiency.

Method	FTG Method	MCS Method	Proposed Method
Computational cost	16,384	1000	680

presents a comparison of the computational efficiency of the MCS and FTG methods and the method proposed in this paper.

In Table 4, the computational cost refers to the number of times the finite element model of the rotor system is invoked during the uncertainty analysis process. For example, the FTG method, which employs a full tensor product mesh with four integration points in each dimension, requires a total of 16,384 integration nodes, necessitating 16,384 invocations of the source model. In contrast, sparse-grid numerical integration requires only 680 invocations. Furthermore, when looking at the percentage reduction in computational cost, sparse-grid integration requires only 68% of the computational cost associated with the MCS. Compared to the FTG method, sparse-grid integration reduces the computational cost by approximately 95.87%, highlighting its superiority in high-dimensional cases. Therefore, employing sparse-grid integration in high-dimensional cases can significantly enhance computational efficiency.

5. Conclusions

In this paper, the method of Chebyshev orthogonal polynomials is utilized to compute the steady-state response of a rotor with uncertainties. To enhance computational efficiency in high-dimensional scenarios involving multiple uncertain parameters, the sparse-grid numerical integration method is employed. A dynamically similar rotor of a turboshaft engine is selected as the study’s subject, and the influences of multiple uncertain parameters on the rotor’s dynamics are analyzed. The main conclusions are as follows:

(1)

The uncertainty in support stiffness significantly affects the steady-state response and critical speed, causing the resonance frequency, originally characterized by a single-peak pattern, to transform into a broad resonance band; this is accompanied by shifts in the resonance peak.

(2)

The upper and lower bounds of the steady-state response are symmetric around the deterministic curve under the influence of mass-unbalance uncertainty, with these bounds widening at high rotational speeds, although the critical speed remains unchanged. Additionally, the phase of the uncertainty unbalance significantly influences the dynamic response of the rotor, resulting in substantial increases in vibration amplitude at critical speeds.

(3)

The effects of multiple uncertainty parameters together on the rotor’s steady-state response are multifaceted, characterized by broad amplitude bands, resonance peaks, frequency shifts, and fluctuations in amplitude magnitude.

In summary, the method proposed in this paper integrates Chebyshev orthogonal polynomial approximation with sparse-grid numerical integration techniques, targeting the analysis of high-dimensional random variables. By utilizing sparse-grid integration, the approach addresses the computational challenges associated with the “curse of dimensionality”, allowing for efficient processing without compromising accuracy. The combination of these techniques not only streamlines the computation but also enhances the robustness of the results, making it a valuable method for engineers and researchers dealing with high-dimensional problems in rotor dynamics and related fields.

Future research will concentrate on the experimental validation of uncertainties in rotor systems to enhance the connections between theoretical models and practical applications. Systematic experimental protocols will be developed to elucidate how various uncertainty factors influence rotor dynamic performance, thereby providing empirical evidence to confirm the accuracy of these models. Furthermore, a strong emphasis on robust optimization methods is crucial, as it integrates results from uncertainty analysis to create designs that are both efficient and resilient to variations in operational conditions.