Robust Dynamic Optimization of a Dual-rotor System Based on Polynomial Chaos Expansion

 Abstract: A novel methodology of robust dynamic optimization of a dual-rotor system based on polynomial chaos expansion (PCE) is developed in this paper. The dual-rotor system model was built by the Timoshenko theory and the finite element method. Instead of the direct Monte Carlo simulation (MCS), the PCE of the present dual-rotor system under support stiffness uncertainty is generated to facilitate a rapid analysis of stochastic responses and yield desirable results in significantly less number of functional evaluations. The PCE is explored as a basis for robust optimization, focusing on the problem of minimizing the unbalance response at operating conditions while minimizing its sensitivity to uncertainty in the support stiffness. This strategy avoids the use of MCS in order to effectively increase the efficiency of the optimization and significantly reduce the computing cost and time spending. The robust dynamic optimization attempts to both optimize the mean performance and minimizes the variance of the performance simultaneously. The multi-objective optimization results show that vibration response can be decreased and is significantly less sensitive to the variation of design parameters compared with initial design case by matching of unbalance amplitude and phase angle differences. Implementation of the proposed robust dynamic optimization in the present dual-rotor system illustrates its potential for further complicated applications.


Introduction
Robust optimization aims to determine the optimal design variables of the system while minimizing the influence of uncertainties [1]. The robust dynamic optimization of rotor systems has attracted the attention of engineers and researchers especially in the fields of the aeronautical and astronautical engineering. Li et al. [2] developed a robust convex optimization for balancing the high-speed rotating machinery. Ritto et al. [3] studied the robust optimization of  Chao-Ping Zang c.zang@nuaa.edu.cn 1 College of Energy and Power Engineering, Nanjing University of Aeronautics and Astronautics, Jiangsu Province Key Laboratory of rotor-bearing systems using Campbell diagram. Penalty functions are introduced to penalize natural frequencies of the systems close to the operational rotational speeds.
Considering the uncertainties of residual unbalances, Stocki et al. [4] developed a robust optimization technique of a typical single-span rotor-shaft using Latin Hypercubus in scatter analysis of the vibration responses. Yang et al. [5] applied the Taguchi method to the dynamic tolerance design of a rotor system and the design of experiment (DOE) is used to make a list of combinations chosen by orthogonal tests. Antinori et al. [6] illustrated a methodology of multidisciplinary robust design process for a fast and flexible low pressure turbine, together with the sensitivity analysis. Hong et al. [7] exploited the response surface method and tolerance model for robust design of a rotor dynamics considering the variation of support stiffness. In recent years, various algorithms such as normal boundary intersection algorithm, the second generation non-dominated sorting genetic algorithm (NSGA-Ⅱ), etc., were used to the robust optimization of rotor system [8][9].
One key component of robust optimization is the quantification of uncertainties in the system output performances propagated from uncertain inputs, named as uncertainty propagation (UP) [10][11]. UP has dominant influence on efficiency and accuracy for robust dynamic optimization of the rotor system. However, as the rotor dynamics design is complicated due to the uncertainties and a significant number of cases based on the traditional Monte Carlo simulation (MCS) are required for evaluating the response statistics of the performance, the robust dynamic optimization process is expensive and time consuming especially when the numerical models are large and complex. Polynomial chaos expansion (PCE) is a newly developed method of UP with attractive attributes for its strong mathematical rigor and ability to produce functional representations of any stochastic quantities. polynomial chaos were to be found in the homogeneous chaos expansion firstly proposed by Wiener [12]. An extension, known as generalized polynomial chaos, was proposed by Xiu and Karniadakis [13]. One important advantage of PCE is that it can conveniently generate the complete probabilistic distribution (e.g., probability density function) of the output response, thus the statistics (mean and variance, etc.) of output performances can be obtained. Recently, the PCE method has been widely studied and applied in engineering, such as model verification and validation [14], aerodynamics [15][16][17] and tolerance prediction [18], etc. In the field of rotor dynamics, Didier et al. [19] investigated the quantification of uncertainty effects on response variability in rotor systems by PCE, but did not perform an optimization. Sinou et al. [20] investigated the influence of PCE order on an uncertain asymmetric rotor system response. Despite the relatively extensive literature available on the application of PCE, few work can be found in the robust design of the rotor system.
In this paper, a robust dynamic optimization method of a dual-rotor system based on PCE is presented. The polynomial chaos expansion is used as a basis for robust optimization, focusing on the problem of minimizing the unbalance response at operating conditions while minimizing its sensitivity to uncertainty in the support stiffness. This strategy avoids the use of Monte Carlo simulation, in which costs can become prohibitive in optimization problems.
The rest of this paper is organized as follows. In Sec. 2, a dual-rotor system is introduced and its modeling is briefly explained for the solution of mechanical problems. In Sec. 3, expansions of the operator of random support stiffness on the chaos basis are explained. The robust dynamic optimization is carried out with the aim to minimize the expected value and variance of unbalanced responses at different operating speeds in Sec. 4,. As a comparison, the optimal solution obtained from the initial design is used to verify the proposed approach. Sec. 5 provides a summary and conclusion.

Dual-rotor System Design Model
The dual-rotor system has been widely used in modern aeroengines and becomes an irresistible technique for increasing the high thrust-weight ratio in aero-engine design process [21]. Figure 1 shows a typical dual-rotor design model which is consisted of high-pressure rotor and low-pressure rotor. The low-pressure rotor adopts 1-1-1 supporting program in which the supports are labeled by number 1, 2 and 5, and the high-pressure rotor adopts 1-0-1 supporting program in which the supports are labeled by number 3 and 4. Because of high rotating speed, rolling bearings are usually adopted and stiffnesses of rotor and support are isotropy. The dualrotor system in the aero-engines has many advantages due to the high and low pressure rotors rotating, reversely. For example, gyroscopic moment, reaction forces of supports and forces conveyed to the aircraft can be greatly reduced. To pursue high thrust-weight ratio, the structure of counterrotating dual-rotor was usually adopted in the rotor system for the advanced aero-engines. where where [M d ] and [G d ] are the mass and gyroscopic matrices of the disc. The vector {f d } defines the unbalance forces due to an eccentric mass. In order to obtain sufficiently reliable results of numerical simulations, together with a reasonable computational efficiency, the rotor subsystem of a rotary machine is usually modeled by means of the one-dimensional finite elements of the beam-type. In the case of only transverse bending vibrations, the rotor is modeled by circular beam elements each having eight degrees of freedom (shown in Figure 2) with rigid discs attached to model compressor and turbine ·4· discs.
The finite element of the Shaft For the degree of freedom [u v α β] T , the vector {f d } is given by where Ue= me×de, me is the mass unbalance and de is the eccentricity. ϕe is the phase shift angle defined by the initial angular position. After assembling the shaft elements and the rigid discs, the equation of motion for the complete rotor system is defined as follows: where Rotor dynamic analysis of the dual-rotor system was performed using the finite element model (FEM) implemented in the in-house code developed in MATLAB. Both the high-pressure rotor and low-pressure rotor are modeled as conservative linearly elastic axisymmetric bodies. The shaft of low pressure rotor is discretized into 22 Timoshenko beam finite elements, and the shaft of high pressure rotor is discretized into 17 Timoshenko beam finite elements. Figure 3 shows a FEM schematic of the dual-rotor system.

Stochastic Model Constructed by PCE under Support Stiffness Uncertainty
The dynamic responses of the actual rotor system are stochastic due to manufacturing inaccuracies related to mechanical properties such as mass, stiffness, damping, and geometric imperfections. In this section, we present a brief introduction to polynomial chaos expansion. Following that we establish the stochastic response model of the present dual-rotor system under support stiffness uncertainty by the PCE method.

The PCE Method
The generalized polynomial chaos expansion of stochastic response Y can be expressed as the following finite number of random variables and a finite expansion order where ξ is the random vector and its dimension is usually consistent with the dimension of the original random input variables in the probability space, ψj(ξ) is an orthogonal polynomial basis. In the Wiener-Askey scheme, different polynomial types associate with different random variables. For more detailed introductions of different pairings of basis polynomials and distributions, readers can refer to Ref. [13]. N is number of the PC used in the expansion excluding the zeroth order term and can be determined by where p is the order of PC and d is the dimension of random vector ξ.
Each of the coefficient αj in PCE (Eq. (5)) can be obtained by the Galerkin projection method for its high accuracy and good robustness. Through projecting both sides of Eq.  (9) where () W ξ is the weighting function corresponding to the polynomial basis ψi(ξ). Taking the expectation of both sides of Eq. (7), αj can be obtained as follow In Eq. (10), the denominator term can be evaluated analytically, while the expectation in the numerator needs to be evaluated with sampling or numerical integration schemes.

Stochastic Model Construction
Considering the influence of support stiffness uncertainty, the elementary stiffness matrix of the support is randomly defined as where kx(τ), ky(τ) are the support stiffness parameters in x and y directions, respectively. kx(τ) is equal to ky(τ) as the support has isotropic property. The random stiffness values of support No. 1, support No.2, support No.3, support No.4 and support No.5 are denoted as k1(τ), k2(τ), k3(τ), k4(τ) and k5(τ) (for convenience, the symbol "τ" will be omitted in the subsequent description). k1, k2, k3, k4 and k5 are mutually independent and their uncertainties are defined in terms of tolerance with 5% (Δki=0.05·ki). It is assumed that k1, k2, k3, k4 and k5 are normally distributed. Using the worst-case formulation, the standard deviations are where σ k i is the standard deviation of random parameter ki.
Then ki can be modelled by the truncated Gaussian random variable ξi, as follow = i i i k i kkσ ξ  (13) where k ̅ i is the nominal stiffness value, ξi is the truncated Gaussian random variable with a zero mean and unit variance. In the probability space Ω(ξ1,ξ2,ξ3,ξ4,ξ5), stochastic response Y can be expanded on the PC basis such that where there is a one-to-one correspondence between ξji and ξi. H(·) is the orthogonal Hermite polynomial basis and H(ξ) of order n is defined from  (15) According to the above, the sequence of processes to establish PCE for the stochastic response under stiffness uncertainty is given in Figure 5.
Random stiffness parameters (k 1 ,k 2 ,k 3 ,k 4 ,k 5 )  Figure 5 The sequence of processes to establish PCE

Stochastic Response
When p=2, Hj(ξ) and E[H j 2 (ξ)] are listed in Columns 3 and 4 of Table 2, respectively, the full factor numerical integration is used to evaluated [Y•H j (ξ)]. The number of quadrature nodes of ξi is set to 3. The quadrature nodes are (-1.7321, 0.0, 1.7321) and the weight values are (0.1667, 0.6667, 0.1667), respectively. The direct tensor product is applied to produce the multi-dimension quadrature nodes and the corresponding weights, resulting in a total of 243 function evaluations to compute the PCE coefficients. Figure 6 shows the stochastic responses at support No.3 (the front fulcrum of the highpressure rotor).  Figure 6 Stochastic responses at support No.3 In Figure 6, the upper and lower envelops (black dashed lines) of the stochastic responses， calculated in a workstation with 32 GB of RAM and two 6-core, Intel Xeon(R) X5690 @ 3.46 GHz processors, are obtained by the PCE method in 238 seconds . The stochastic responses (gray solid lines) are obtained by the MCS, in which 1000 samples are chosen. The CPU time is about 1126 seconds, which is around 4 times higher than that of PCE method.

Robust Dynamic Optimization
Robust dynamic optimization is inherently a multi-object optimization problem in which both the expected value(s) and the standard deviation(s) of the performance function(s) are to be minimized [22]. The different types of factors and their relationships with the rotor system can be represented by P-diagram, as shown in Figure 7. The objective is to minimize the unbalance response and variance at each operating condition. For the present dualrotor system, the operating speeds ωi=1000+i·1000 r/min (i =1,2,…,8) of the high-pressure rotor are considered. The first objective function is where s={ω1, ω2, …, ω8} is the operating speeds. x={x1, x2, x3} is the design parameters. x1 and x2 denote the unbalance masses placed at HPC4 and HPC9; x3 denotes the phase shift angle of HPC9. Δk={Δk1, Δk2, Δk3, Δk4, Δk5} represents the uncertainty of support stiffnesses. μ v i (·) is the mean of unbalance response of support No.3 at operating speed ωi (i=1,2,…,8). The second objective function is where σ v i (·) is the standard deviation of unbalance response of support No.3 at operating speed ωi. Different ranges of the design variables are shown in Table 3.  Figure 8 shows the Pareto fronts (Pareto front-a set of solutions such that no objective can be improved further without worsening another at the same time) for all cases shown in Table 3.

Figure 8 Pareto fronts
As shown in Figure 8, the mean and variance of the response at support No.3 tend to decrease as the low boundary becomes smaller. However, the smaller unbalance mass is, the higher the requirement for machining and balancing precision is. Different decisions can be specified based on realistic goals. For illustration, the Pareto front for Case 6 (the low boundary equals to 50 g·mm) is shown in Figure 9. Table 4 lists the initial design and selected robust design RD1, RD2 and RD3 in Case 6.  From Figure 9, RD3 may be selected when the variance is Ben-sheng Xu et al.

·8·
expected to be as small as possible. In Table 4, InitD represents the initial design. It can be seen that the vibration response of the support does not decrease simply with the decrease of the amplitude of unbalance on HPC9. Numerical calculations show that a significant reduction of the vibration level and robust design under support stiffness uncertainty is obtained if phase angle is considered as design parameter. To demonstrate the effective of the robust design approach, MCS is used to evaluate the possible response variations of the rotor system under the stiffness uncertainty. Figure 10 shows mass unbalance responses at support No.3 for InitD and RD1. In Figure 11, RD1 (green solid lines) and InitD (gray solid lines) are investigated and frequency response variations at support No.3. It can be seen that the possible maximum response (about 82.9 mm/s) under the stiffness uncertainty in RD1 is over ten lower than that of InitD (the possible maximum response is about 955.0 mm/s).

Conclusions
In this paper, we have described a novel robust dynamic optimization method of a dual-rotor system based on PCE.
The dual-rotor system model was built by the Timoshenko theory and the finite element method. In the presence of support stiffness uncertainties that is defined in terms of tolerance, 243 function evaluations were performed to compute the PCE coefficients. The PCE was generated to facilitate a rapid analysis of stochastic responses and the results are similar to what were obtained using 1000 Monte Carlo samples with more efficient calculations. As a result, the PCE bases utilized make the present dual-rotor design model generalized for stochastic computations and reduce computational effort. The robust dynamic optimization (RDO) of the dual-rotor system was performed by minimizing the expected value and variance of unbalanced response at different operating speeds. The efficiency of the proposed RDO method was tested and validated through Monte Carlo simulation of the effects of uncertainties. The multi-objective optimization result is shown to lead to designs that significantly decrease vibration responses of the rotor system and improved system performance with reduced sensitivity to the given uncertain stiffness parameters. The study has paved way for the proposed PCE based RDO method to be applied to further complex large-scale applications.

Declaration
Funding Supported in part by the National Natural Science Foundation of China and National Safety Academic Foundation of China (Project No. 12072146, U1730129), the supports from Jiangsu Province Key Laboratory of Aerospace Power System, the Key Laboratory of Aeroengine Thermal Environment and Structure, Ministry of Industry and Information Technology are also gratefully acknowledged.

Availability of data and materials
The datasets supporting the conclusions of this article are included within the article.

Authors' contributions
The conceptualization and methodology were developed by Xu and Zang. The original draft of the manuscript was written by Xu and the review and editing was made by Zang. Wang and Zhang joint a constructive discussion of the content. Zang is in charge of corresponding the paper.