A Reduced-Order Model for the Vibration Analysis of Mistuned Blade–Disc–Shaft Assembly

: An e ﬀ ective reduced-order model is presented in this paper for the vibration analysis of a mistuned blade–disc–shaft assembly considering the ﬂexibility of the shaft and the rotordynamic e ﬀ ects. For the sake of accurate modeling and quantitative analysis, three-dimensional (3D) ﬁnite element models were employed in obtaining the governing equations of motion with the Coriolis force, centrifugal sti ﬀ ening, and spin softening e ﬀ ects taken into account. Then, an e ﬃ cient model order reduction technique based on the coordinate projection by normal modes of tuned assembly and cyclic symmetry analysis was developed for mistuned blade–disc–shaft assembly. The criterion of whether one matrix could be incorporated in cyclic symmetry analysis is presented. During the modeling, the mistuning in blade and disc was taken into account and dealt with independently. In mistuning projection, the blade and disc parts were both projected onto their tuned counterparts of the sector model, where the boundary conditions were set to be ﬁxed and free, respectively. Finally, an example of a blade–disc–shaft assembly was employed to validate the e ﬀ ectiveness of the presented method in free and forced vibration analysis.


Introduction
Blade-disc-shaft assemblies are the key components of engines and compressors and take an important position in the design of these machines [1]. As the basis for design, developing accurate dynamic models for such assemblies and gaining deep insight into their vibration characteristics are of great significance. For this reason, the vibration issues of blades and rotor systems attracted enormous attention from researchers and engineers, and a huge number of investigations were published in the past few decades [1][2][3][4]. Nevertheless, most of the available investigations concentrated either on the vibration of the blade and disc by neglecting the flexibility of the shaft or on the dynamics of rotor systems by assuming the blade and disc to be rigid and tuned.
Some researchers considered the mutual interactions between bladed discs and shafts and developed several dynamic models for blade-disc-shaft assemblies. Crawley et al. [5] presented a simplified analytical model for rotating a flexible blade-rigid disc-flexible cantilevered shaft system. Okabe et al. [6] addressed an equivalent modeling technique for a blade-shaft system by dividing the system into a shaft plus additional equivalent mass-spring models of blades. Zhang et al. [7] presented a method for dynamic analysis of flexible blade-disc-shaft systems by means of the finite element method and a modal synthesis approach. Huang et al. [8] developed an analytical approach for the coupled shaft torsion and blade bending vibration analysis using the weighted residuals method. Chun et al. [9] proposed an analytical substructure synthesis method for analyzing the coupled free vibrations between shaft and blades. A1-Bedoor [10,11] presented efficient models for coupled shaft-torsional and blade-bending vibration analysis using Lagrange's equation and the assumed modes method. Ma a 3D finite element model of a mistuned blade-disc-shaft assembly, its governing equation of motion at a constant rotating angular velocity Ω has the following general form: where M, C d , and K e are the mass, damping, and elastic stiffness matrices of the assembly, C c (Ω), K c (Ω), and K s (Ω) represent the speed-dependent Coriolis, centrifugal stiffening, and spin softening matrices, respectively, K δ denotes the additional stiffness matrix due to mistuning, u is the vector of response, and f represents the vector of external forces. As the 3D finite element models are employed, the above matrices are obtained via linear assumption. In this paper, the external forces acting on the blades are assumed to be traveling wave excitations, where the forces on the ith blade can be expressed as where B 0 denotes the forcing amplitude vector, ω represents the angular frequency of excitation, and ω = CΩ, where C is the engine order (EO) of excitation, and ϕ i is the interblade phase angle defined as ϕ i = (i − 1)Cθ 0 , where θ 0 = 2π/N, and N is the number of blades.

Model Order Reduction via Cyclic Symmetry Analysis
Model order reduction techniques are often employed, combined with finite element models, to significantly reduce the computational requirements during analysis. It is much more urgent for mistuned bladed discs, as mistuning is highly random and thousands of mistuning patterns are often needed to study the statistical vibration characteristics. In this paper, the ideas presented in References [28][29][30] are employed and further developed for blade-disc-shaft assemblies. From a general perspective of view, these methods belong to the mode superposition method, whereas the modes used are those of a tuned bladed disc rather than a certain mistuned case. According to these approaches, the normal modes of the tuned blade-disc-shaft assemblies are used as the basis for projecting the mistuned physical model into a modal domain, which can be represented as where Φ t denotes the mode matrix of the tuned assembly in static condition, where the matrices due to rotation are not included during eigenvalue analysis, and q represents the modal coordinate vector. The effectiveness of this coordinate transformation for mistuned blade-disc-shaft assemblies is validated in the next section.
In the vibration analysis of bladed discs, the cyclic symmetry properties are often employed to reduce the computation requirements [31,32]. A tuned blade-disc-shaft assembly is also cyclically symmetric if the shaft is axially symmetric. In the analysis of a rotating blade-disc-shaft assembly, the rotordynamic effects should not be neglected. However, few relative investigations considered these effects. As the rotordynamic matrices are all speed-dependent, the Coriolis matrix is skew-symmetric, and the centrifugal stiffening matrix depends on the stress distributions due to centrifugal forces, whether these matrices can be included in cyclic symmetry analysis should be rigorously verified. It can be proven that one matrix could be incorporated in cyclic symmetry analysis if such a matrix of the sector model has the following relationship with that of its adjacent sector. where M s and M s+1 denote the matrices of sth and (s + 1)th sectors, respectively, and R s (θ 0 ) represents a rotation transformation matrix in term of the sector angle.
where I N n denotes an identity matrix of dimensions N n , where N n is the number of nodes in a sector model, the symbol ⊗ denotes the Kronecker product, and R e (θ 0 ) is an elementary matrix. Taking the rotation axial along the z-axis as an example, R e (θ 0 ) can be expressed as Detailed proof of this proposition, for brevity, is presented in Appendix A. According to the explanation in Appendix A, the Coriolis, centrifugal stiffening, and spin softening matrices can all be incorporated in cyclic symmetry analysis. In this circumstance, the mode matrix Φ t represented in terms of cyclically symmetric modes can be used as the projecting basis during model order reduction. Then, Φ t can be explicitly expressed as where E denotes the complex Fourier matrix with e mn = e j(m−1)(n−1)θ 0 / √ N for m, n = 1, 2, · · · , N, I N s denotes an identity matrix, Φ s h is the mode matrix of hth nodal diameter (ND), and Bdiag where where Appl. Sci. 2019, 9, 4762 5 of 19 where the subscripts of "1", "2" and "3" in the above matrices represent the terms in the corresponding block-circulated matrices in Equation (A10) shown in Appendix A. The reduced-order mistuning matrix K δ is also generated via the coordinate transformation in Equation (3), giving In most available investigations, only the mistuning in blades was taken into account. This may partly have been because the blade tends to be more flexible and less structurally coupled than the disc. Actually, mistuning in discs is also inevitable due to manufacturing tolerances and deteriorations during operation. The mistuning in discs may not play that significant a role as blades in the vibration of bladed discs. Nevertheless, it may be non-negligible for blade-disc-shaft assemblies, as the flexibilities in the shaft may lead to significant coupling motions between the shaft and disc. Therefore, the mistuning in blades and discs is considered and dealt with independently in this paper. Then, K δ can be further expanded, without loss in accuracy, into blade and disc parts as where Φ b t and Φ d t denote the mode matrices corresponding to the blade and the disc parts, respectively, and K b δ and K d δ represent the corresponding mistuning matrices of the two parts. Despite the simple expressions, K b δ and K d δ are difficult to determine, as the mistuning matrices K b δ and K d δ of a realistic bladed disc are highly random and difficult to measure. For this reason, frequency mistuning models were more often employed in available investigations for mistuning quantization. In this paper, such a model is also used. Detailed deductions on obtaining the reduced-order mistuning matrices K b δ and K d δ are presented in the sections below.

Parametric Projection of Blade Mistuning
The mistuning projection method proposed by Lim et al. [30] is employed in this section to include blade mistuning in the reduced-order model. The core idea of this method is to project the mistuning matrix onto the normal modes of a tuned blade cantilevered at its root. The main procedures of this method are briefly reviewed herein. Firstly, the blade part mode Φ b t is represented by cyclically symmetric modes as shown in Equation (7).
where Φ Then, the blade part's cyclically symmetric modes are represented by the normal modes of the tuned blade cantilevered at its root.
where Φ b f ix is the normal mode matrix of the tuned cantilevered blade, and p b h denotes the participation factor. By the above transformation, Appl. Sci. 2019, 9, 4762 δ,h tends not to be diagonal, it is often approximately processed in frequency mistuning models by keeping only the diagonal elements, which can be represented as where M is the number of retained modes of the tuned cantilevered blade, λ b m represents the mth eigenvalue of the tuned cantilevered blade, and δ b h,m denotes the deviation ratio, which is often set to be a random parameter during analysis.

Parametric Projection of Disc Mistuning
Unlike the projection of blade mistuning, the coordinates between adjacent sectors of the disc are directly coupled, and the disc mistuning matrix cannot be expressed in block-diagonal form. Therefore, the mistuning projection method in the previous section cannot be directly employed. Nevertheless, the idea of representing the overall modes by normal modes of a relatively simple component and projecting the mistuning matrix onto these normal modes can be utilized. Of course, some modifications to the formulations are needed to account for the coupling issue of the disc model.
Firstly, the disc part mode Φ d t is represented by cyclically symmetric modes.
where Φ d h denotes the disc part of Φ s h , and I N d is an identity matrix of dimensions N d , where N d is the number of DOFs in a sector model of the disc with the left cyclic symmetry surface eliminated.
As mentioned above, the mistuning matrices of adjacent sectors of the disc are crossed and coupled. Proper processes are needed to make the matrix block-diagonal. It can be found that adjacent sector models can be made independent in form by coordinate expansion, where the original and expanded coordinates with respect to the disc model are expressed as where u d t denotes the original coordinate vector with respect to Φ d t ,û d t is the expanded coordinate, and the subscript "R", "L" and "I" denote the right and left cyclic boundary DOFs and internal DOFs, respectively. It can also be proven that the disc mistuning matrix K d δ can be obtained using an expanded block-diagonal matrixK where where R d (θ 0 ) denotes a rotation transformation matrix, and R d (θ 0 ) = I cd ⊗ R e (θ 0 ), where I cd is an identity matrix of dimensions N cd , and N cd is the number of nodes on cyclic symmetry interface.
Then, similar to that shown in Equation (14), the expanded cyclic modesˆ Φ where Φ d f re denotes the normal mode matrix of a tuned sector model of the disc with the left and right interfaces, as well as the adjacent interface between the disc and shaft, being set to be free, and p d h denotes the participation factor. Herein, the free constraint is employed on the disc-shaft adjacent interface instead of a fixed one, because the latter one may not accurately describe the 1ND coupling modes between shaft and disc, leading to non-negligible errors in forced responses. The employment of free modes avoids such issues and presents satisfactory results in both frequencies and forced responses.
By substituting Equation (22) into the expression of K d δ and neglecting the non-diagonal terms, the disc part's reduced-order mistuning matrices K d δ can be expressed as where the deviation ratio, and λ d w denotes the wth eigenvalue of the tuned disc component. On this basis, the free and forced vibration of tuned and mistuned blade-disc-shaft assemblies can be analyzed by using the governing equation of motion in Equation (8) and the coordinate transformation in Equation (3). As the mistuning in blade and disc is taken into account, the effects of blade and disc mistuning on the vibration characteristics of blade-disc-shaft assemblies can be investigated.

Validation of Reduced-Order Models
Validation of the developed order reduction method is presented in this section. As the accuracy of the reduced-order model in vibration analysis for practical structures relies on the accuracy of its counterpart full-order model, the effectiveness of the developed method can be verified by comparing the results of the reduced-order model with those of the full order model. Figure 1 illustrates the finite element model of a blade-disc-shaft assembly, which is used as an example in this section, where the DOFs of the whole model and the sector model are equal to 103,563 and 7179, respectively. The whole model contains 18 sectors. The reduced-order model contains 216 DOFs with 12 modes per ND. During the projection of blade and disc mistuning, the retained modes of the cantilevered blade and free disc are 20 and 50, respectively. The excitation is applied at the tip of the blades with an amplitude of 2 N. The structural damping is assumed to be Rayleigh damping and C d = βK e , where β is set to 3 × 10 −6 in the subsequent analysis. The material of the blade and disc is titanium alloy, and that of the shaft is steel, where the properties of these materials are listed in Table 1. retained modes of the cantilevered blade and free disc are 20 and 50, respectively. The excitation is applied at the tip of the blades with an amplitude of 2 N. The structural damping is assumed to be Rayleigh damping and = d e  C K , where  is set to 3 × 10 −6 in the subsequent analysis. The material of the blade and disc is titanium alloy, and that of the shaft is steel, where the properties of these materials are listed in Table 1.  Generally, the mistuning in the blade and disc is random and hard to measure. Herein, a simpler mistuning model, the commonly used proportional mistuning model, is used to validate the developed method. Small random variations in Young's modulus are introduced as  Table 2.   Generally, the mistuning in the blade and disc is random and hard to measure. Herein, a simpler mistuning model, the commonly used proportional mistuning model, is used to validate the developed method. Small random variations in Young's modulus are introduced as where E b n and E d n denote the Young's modulus of the nth sector models of the blade and disc, E 0 represents the nominal Young's modulus, and δ b n and δ d n are random variables representing the small variations. Specifically, the employed proportional mistuning patterns of the blade and disc in the subsequent analysis are listed in Table 2.  Figure 2 illustrates the natural frequencies versus nodal diameters of the tuned blade-disc-shaft assembly in a static condition. Some of the 0ND and 1ND modes are shaft-dominated and not involved in the lines. With this diagram, the excited modes can be easily obtained by comparing the exciting frequency and the natural frequencies.
Appl. Sci. 2019 Figure 2 illustrates the natural frequencies versus nodal diameters of the tuned blade-discshaft assembly in a static condition. Some of the 0ND and 1ND modes are shaft-dominated and not involved in the lines. With this diagram, the excited modes can be easily obtained by comparing the exciting frequency and the natural frequencies.  Figure 3 shows the comparisons of the forced responses of the mistuned blade-disc-shaft assembly using the normal modes of the mistuned case and those of its tuned counterpart. The presented results are the responses of the excitation nodes, which are located at the tip of each blade. In the subsequent analysis, the results of these nodes are also employed. It can be seen in the figures that the results agree very well. The maximum relative error in amplitude is less than 0.1%. Therefore, the normal modes of the tuned blade-disc-shaft assembly can accurately represent the modes of mistuned assembly.   Figure 3 shows the comparisons of the forced responses of the mistuned blade-disc-shaft assembly using the normal modes of the mistuned case and those of its tuned counterpart. The presented results are the responses of the excitation nodes, which are located at the tip of each blade. In the subsequent analysis, the results of these nodes are also employed. It can be seen in the figures that the results agree very well. The maximum relative error in amplitude is less than 0.1%. Therefore, the normal modes of the tuned blade-disc-shaft assembly can accurately represent the modes of mistuned assembly.  Figure 4 shows the comparisons of the first 100 natural frequencies of the tuned assembly using the full-order model and cyclic symmetric sector model. It can be seen in the figure that the two cases agree well, with the maximum relative error in natural frequency being less than 0.01%. Therefore, cyclic symmetry analysis is also feasible for bladed discs mounted on flexible shafts.    Figure 4 shows the comparisons of the first 100 natural frequencies of the tuned assembly using the full-order model and cyclic symmetric sector model. It can be seen in the figure that the two cases agree well, with the maximum relative error in natural frequency being less than 0.01%. Therefore, cyclic symmetry analysis is also feasible for bladed discs mounted on flexible shafts.  Figure 4 shows the comparisons of the first 100 natural frequencies of the tuned assembly using the full-order model and cyclic symmetric sector model. It can be seen in the figure that the two cases agree well, with the maximum relative error in natural frequency being less than 0.01%. Therefore, cyclic symmetry analysis is also feasible for bladed discs mounted on flexible shafts.    Figure 5 shows the relative errors in eigenvalues of the reduced-order model at different rotating speeds, where the first 20 are included and the mistuning in the blade and disc is considered. It can be seen in the figure that the relative error grows with rotating speed but within a relatively small range. Although the relative errors increase for higher-order eigenvalues, the overall maximum relative error is still less than 0.1% when the rotating speed is less than 2000 rad/s. Thus, the presented model order reduction technique can give accurate results in eigenvalues for a mistuned blade-disc-shaft assembly in spite of the significant reduction in model order. Figure 5 shows the relative errors in eigenvalues of the reduced-order model at different rotating speeds, where the first 20 are included and the mistuning in the blade and disc is considered. It can be seen in the figure that the relative error grows with rotating speed but within a relatively small range. Although the relative errors increase for higher-order eigenvalues, the overall maximum relative error is still less than 0.1% when the rotating speed is less than 2000 rad/s. Thus, the presented model order reduction technique can give accurate results in eigenvalues for a mistuned blade-disc-shaft assembly in spite of the significant reduction in model order.    In order to illustrate the efficiency of the presented method, the computational time of the reduced-order model in computing natural frequencies and forced responses is compared with that of the full-order model, which is listed in Table 3. The specifications of the employed computer during the analysis were as follows: Intel Core i7-8700 central processing unit (CPU) with 3.20 GHz; size and speed of the cache memory = 32 GB and 2666 MHz, respectively; size of the solid-state drive (SSD) hard disk = 512 GB. The presented values are all the averages of 100 analysis cases. It can be seen that the presented model order reduction technique could significantly increase the computational efficiency. Table 3. Comparison on the computational efficiency between full-order model (FOM) and reduced-order model (ROM). Based on the results presented in this section, the effectiveness and efficiency of the presented method can be validated. With this method, the coupling vibrations between bladed disc and shaft, as well as the effects of the shaft's flexibility on the vibration of a mistuned bladed disc, can be studied. Nevertheless, these issues are not the focus of this paper and can be addressed in further studies.

Conclusions
In this paper, an effective reduced-order model was presented for the vibration analysis of a mistuned blade-disc-shaft assembly. Three-dimensional finite element models, which can provide accurate modeling and quantitative analysis for complex structures, were introduced to obtain the governing equations of motion, where the Coriolis force, centrifugal stiffening, and spin softening effects due to rotation were all taken into account. An efficient model order reduction technique, based on the coordinate projection by normal modes of tuned assembly and cyclic symmetry analysis, was developed for a mistuned blade-disc-shaft assembly. As the basis for employing cyclic symmetry analysis, the criterion of whether one matrix could be incorporated in the analysis, as well as rigorous proof, was given. It was shown that the Coriolis, centrifugal stiffening, and spin softening matrices satisfy this criterion and can be converted into block-circulated forms in cyclically symmetric coordinates. In the presented model order reduction method, the mistuning in the blade and disc was taken into account and dealt with independently to reflect the realistic conditions of blade-disc-shaft assemblies. Regarding the coupling issue of the disc mistuning matrix, an effective mistuning projection approach was developed by using coordinate expansion and projecting the mistuning matrix onto a free tuned sector model of the disc. Detailed deductions on the theoretical basis of coordinate expansion were presented. The employed example of a blade-disc-shaft assembly illustrates the effectiveness of the presented method in free and forced vibration analysis in both static and rotating conditions.
As cyclic symmetry analysis was employed in this paper, the presented method required the rotating systems to be cyclically symmetric. In some circumstances, such requirements may not be satisfied due to the anisotropy in supports, such as the anisotropy in bearings and foundations, which may lead to periodically time-variant terms and more complicated vibration characteristics. Then, further development of this method should be done, which was not included in this paper and remains to be conducted in the future.

Conflicts of Interest:
The authors declare no conflicts of interest.

Appendix A
In this section, detailed proof of the proposed proposition in Section 3.1 is presented. Regarding the issue of including the rotordynamic matrices in cyclic symmetry analysis, the determination factors are whether the matrices can be transferred into block-circulated form in cyclically symmetric coordinates. Thus, the core content of the proposition is that one matrix can be transferred into block-circulated form in a cyclically symmetric coordinate system if it satisfies the relationship shown in Equation (4). As an example, a block-circulated matrix is shown herein as follows: where Z 1 , Z 2 , and Z 3 are the sub-blocks of Z. Without loss of generality, an N-sector finite element model is used as an example for detailed deduction. For brevity, the DOFs on the right and left cyclic boundaries are represented by "R" and "L" and the others as "I", and the left part of the sth sector u (s) L is equal to the right part of the (s + 1)th sector u (s+1) R . Then, the coordinates of this model in a global Cartesian coordinate system can be represented as (1) u (2) . . . u (s) . . .
Then, a certain structural matrix Y of this model in a global Cartesian coordinate system has the following form: For clarity, the above matrix is further expressed as block form as 2 , and Y (s) 3 for s = 1, 2, · · · , N denote 2 × 2 block matrices, and The cyclically symmetric coordinate system is equivalent with several independent local Cartesian coordinate systems attached to each sector. Thus, the coordinate transformation between the global Cartesian coordinate system and the cyclically symmetric coordinate system can be expressed as (1) u (2) u (3) u (4) . . . (1) u (2) . . . u (n) . . .
where R(sθ 0 ) for s = 1, 2, · · · , N − 1 represents the rotation transformation matrix in terms of the angle sθ 0 , and I S is an identity matrix. It can be easily understood that the rotation transformation matrix R(sθ 0 ) is equal to the products of R(θ 0 ) s times, i.e., R(sθ 0 ) = [R(θ 0 )] s . For brevity in expression, R(θ 0 ) is simply represented by R 0 , and (R 0 ) N = I S . By applying this coordinate transformation to Y, such a matrix in a cyclically symmetric coordinate system becomes Y = R T c YR c , and If such a matrix satisfies the relationship in Equation (4), the following equation can be obtained for s = 1, 2, · · · , N − 1: With this equation, it can be inferred that Thus, the matrix Y can be expressed as a block-circulated form as follows: where 3 . According to the physical meaning of Coriolis force, the Coriolis and spin softening matrices satisfy the proposition in Equation (4), as these two matrices can be directly generated by finite element models. On the other hand, the centrifugal stiffening matrix depends on the stress distributions due to centrifugal forces. Fortunately, the centrifugal force is a 0ND excitation, which leads to identical deformations and stresses in the local cylindrical coordinate systems of each sector. As a result, the centrifugal stiffening matrices of each sector are also identical in local coordinate systems. Then, the overall centrifugal stiffening matrix possesses also the block-circulated form in a cyclically symmetric coordinate system. Therefore, the Coriolis, centrifugal stiffening, and spin softening matrices can all be expressed as block-circulated forms and be included in cyclic symmetry analysis. (A11) Due to the coupling between sectors, the mistuning matrix K d δ contains crossed terms, whereK d δ has an block-diagonal form.