An Efficient Dynamic Modeling Technique for a Central Tie Rod Rotor

Compared with the three-dimensional rotor model for a central tie rod rotor, an equivalent one-dimensional model can greatly improve the computational efficiency in rotor dynamics analysis with a certain accuracy. However, little research work can be found on improving the modeling accuracy of one-dimensional models using experimental data. In this paper, a onedimensional discrete mass model considering pretightening force is proposed for central tie rod rotors to achieve the purpose of both efficient and accurate modeling. Experimental testing and three-dimensional model analysis are used as reference and verification approaches. A sensitivity-based method is adopted to update the proposed one-dimensional model via minimizing the error in the critical speed comparing with the corresponding three-dimensional finite element model which has been verified by a modal test. Prediction of damped unbalanced response is conducted to show the practicality of the updated onedimensional model. Results show that the method presented in this research work can be used to simulate a complex preloaded rotor system with high efficiency and accuracy.


Introduction
A central tie rod rotor is a commonly used structural form for a gas generator rotor in a turboshaft engine. It is a kind of combined rotor, which has lighter weight and more convenient to manufacture than integrated rotors. The structure of a real gas generator rotor is complex which contains numerous connection structures and assembly structures. The cross-section of a typical central tie rod rotor is shown in Figure 1. The rotor consists of a three-stage compressor disk, a one-stage centrifugal impeller disk, a one-stage gas turbine disk, and a central tie rod. All the disks are compressed axially by a central tie rod through a compression nut. Research indicates that the influence of pretightening force should be considered when the dynamic characteristics of the central tie rod rotor are studied [1].
Three-dimensional (3D) finite element models are often adopted in the rotor dynamics analysis for gas generator rotors [2]. The key factor in modeling a rotor is to accurately simulate its mass, inertia moment, and bending stiffness properties [3]. However, the 3D model is extremely timeconsuming. Therefore, the 3D rotor model is commonly simplified into a one-dimensional (1D) model for analysis. Metsebo et al. [4] investigated the dynamics of a rotor-bearing system using the Timoshenko beam. Nonlinearity of ball bearings was studied through frequency response and bifurcation. In order to improve the computational efficiency and retain the characteristics of a 3D model, Carrera and Zappino [5] proposed a Node-Dependent Kinematic method in which different beam theories can be applied to different nodes in the same beam element. Filippi et al. [6] extended the Node-Dependent Kinematic method based on 1D element; the dynamic characteristics of a cantilever rotor and a multidisk gas turbine rotor were studied, respectively. The research works on the rotor dynamics described above were all based on reduced dimension models. However, the accuracy of these simplified models cannot be guaranteed due to the lack of experimental data.
A key issue in the dynamic analysis of a central tie rod is the influence of pretightening force. Many models of tied-rod rotors considering pretightening force have been developed. Gao et al. [7] used beam elements to simulate the tie rod. The positive pressure in the rod was obtained under pretightening force, and then, the equivalent bending stiffness of the rotor was deduced. The effect of pretension on bending stiffness was well verified by the test results. Meng et al. [8] obtained the equivalent stiffness of the contact area on the disk from the finite element model, and then, the dynamic characteristics of a central tie rod rotor were obtained via transfer method. The effect of pretightening force on the equivalent stiffness of contact areas was investigated. Shaposhnikov and Gao [9] established both the 3D model and the 1D beam model for a central tie rod rotor via finite element software. The error between the two models was reduced via the modification on the stiffness of the contact region according to the strain energy distribution in each mode.
The mathematical model of the structure is always different from the actual system due to assumptions on simplification [10]. Consequently, it is necessary to determine the mechanical parameters to obtain an accurate mathematical model for the simulation [11,12]. Model updating is a technique for determining structural mechanical parameters by using modal frequencies, mode shapes, and frequency response functions [13,14]. Sensitivity analysis is usually used to select model correction parameters [15]. In order to improve the prediction accuracy of the model, model updating has been widely used in engineering, especially in aeroengine modeling [16]. Research works on model updating of rotor systems have been done in recent years. Ricci et al. [17] carried out the model updating of a steam turbo generator. Moment of inertia and stiffness coefficient in the model were modified according to the minimization of the errors in natural frequencies. Chouksey et al. [18] modified the support stiffness and damping of a single disk rotor using the inverse eigensensitivity method, and these parameters were verified by unbalanced response. Later, Chouksey et al. [19] updated the model of the rotor by using the eccentricity of journal bearing and material damping coefficient of the shaft. Thus, the correction parameters of the bearing are reduced from 8 to 1, and a more accurate result is obtained. Jha et al. [20] modified the support stiffness and damping of the Timoshenko beam rotor by a nonlinear least-square optimization method. Frequency response function and unbalanced response were adopted to update the finite element model of a rotor.
Finite element modeling methods were adopted in the aforementioned works. Another kind of discrete model which regards as the lumped mass model is also commonly used. The lumped mass model simplifies the actual rotor into a series of massless elastic shafts and a rigid disk with mass consolidated on the shafts. Compared with the continuous mass model and the finite element model, the lumped mass model has simpler structure that requires less calculation. When reasonable simplification is conducted, the lumped mass model can be accurate in response prediction and much more efficient, which provides a power tool for performing the rotor dynamics analysis on complex rotating machinery.
Although the beam model has been widely adopted in the design stage of a rotor, this model is more rigid than the real rotor [9]. According to the author's knowledge, limited studies on the model updating of the 1D model of the central tie rod rotor based on damped critical speed have been conducted. Therefore, a central tie rod rotor model with pretightening force is established in this paper using the discrete mass model. The 1D rotor model is modified by a sensitivity-based model updating technique to achieve an efficient and accurate prediction of the dynamic characteristics.
Research works in this paper are organized as follows. The simplified modeling of a central tie rod rotor considering the pretightening force is proposed at first in Section 2. Then, a brief introduction on the theory of model updating technique is presented. A free-free modal test of a real central tie rod rotor will be conducted to verify the fine-meshed 3D model referenced for the subsequent model updating. Then, 2 International Journal of Aerospace Engineering a sensitivity-based model updating on the proposed model will be conducted in Section 3. The accuracy of the updated model will be verified by the results of unbalanced response. Conclusions will be drawn in Section 4. The implementation procedure of the methodology proposed in this paper is summarized in Figure 2. The whole procedure mainly includes four steps: (1) simplification, (2) modeling, (3) modification, and (4) response calculation. Data for the validation of the updated model comes from commercial finite element software ANSYS.

Structural Simplification of a Central Tie Rod Rotor.
When two adjacent disks are compressed only due to the axial force, the axial position will retain the same position before and after compression; therefore, it will be regarded as fixed. The drum between two disks is simplified as a cylinder; the disks at all levels are regarded as thin disks without considering the thickness; since aerodynamic effects are not taken into account here, the blades are deleted; the drum shaft between disks is simplified as an elastic shaft regardless of mass; the central tie rod is simplified as an elastic shaft considering mass, and the elastic support is simplified as a spring. According to the above simplification principle, the mass distributions of disc drum and central tie rod are shown    The disk drum and central tie rod share node 1 and node 12. k a and k b represent the support stiffness of the system.

Governing Equations.
In this section, the Lagrange equation is used to establish the differential equation of motion of the rotor. The stiffness of the elastic supports of the rotor is k a and k b , respectively. The distance between the i th disk and the left fulcrum is a i . The distance between the i th mass point on the central tie rod and the left fulcrum is b i . The span between two elastic supports is l. The mass, polar moment of inertia, and equatorial moment of inertia of each point are m i , J pi , and J di , respectively. Each lumped mass point contains four degrees of freedom. x i and y i are the translations along the x-and y-direction at the i th point, respectively. θ xi and θ yi are the rotations along the x-and y-direction at the i th point, respectively. Axial displacement is not considered. The rotation angular velocity of the rotor is denoted by Ω.
The kinetic energy of the rotor includes the rotational kinetic energy and the translational kinetic energy at each center of mass. The kinetic energy of the rotor T r can be expressed as The potential energy of the rotor V r can be expressed as where Symbol K in Equation (4) represents the stiffness matrix of the whole rotor system in a plane considering the influence of elastic support and the pretightening force in central tie rod. It is derived from the inversion of the flexibility matrix. Since the drum rotor is simplified to a stepped hollow shaft, the bending deformation of each drum rotor is calculated by a continuous piecewise independent integral method.
For a compressed segment as shown in Figure 4, assume that the microsegment of the bar is subjected to end forces and external load. According to the shear equilibrium equation, we have where F s ðxÞ is the shear force on the cross-section. qðxÞ is the distribution force per unit length.
If qðxÞ equals to zero, then where MðxÞ is the moment on the cross-section. F N is the axial force. wðxÞ is the deflection of the rod end. From Equation (7), the following equation can be obtained:  Eliminating high order small quantity, taking the second derivative of x, we have For straight beams with equal cross-section, the relationship between bending moment and deformation is where E and I are the elastic modulus and the moment of inertia of the section, respectively. Substituting Equations (6) and (10) into Equation (9), we have The general solution of Equation (11) is For the stretched segment, changing the F N in Equation (7) into -F N , the deductions are the same as those for the compressed microsegment. The general solution of the deflection for the stretched segment is as follows: The coefficients in both Equation (12) and Equation (13) are obtained via force boundary conditions and displacement boundary conditions between stepped shaft sections. Then, the stiffness matrix K can be derived from the inverse of the flexibility matrix.
Substituting kinetic energy and potential energy into the Lagrange equation: By taking the generalized force as zero, the steady-state eddy differential equation of the rotor can be obtained: where m n is the mass of the n th mass point, and J dn is the diameter moment of inertia of the n th mass point.
where J pn is the polar moment of inertia at the n th mass point.
The displacement in the xoz plane can be expressed as

International Journal of Aerospace Engineering
The displacement in the yoz plane can be expressed as The system damping is assumed to be proportional damping, so the final governing equation can be obtained.
where M r is the mass matrix, G r is the gyroscopic matrix, K r is the stiffness matrix, and F r is the generalized force vector. The damping of the system C r is considered as proportional damping, C r = αK r + βM r , and the damping coefficient α and β can be measured experimentally. The algorithm and simulation proposed in this paper are based on MATLAB. Thus, damped critical speed and unbalanced response can be obtained.

Model Updating Based on Sensitivity.
The simplification on the reference model can greatly improve the efficiency during the rotor dynamics analysis, yet it will lead to errors in the analyzed results, so it is necessary to modify the design parameters of the rotor to achieve a good accuracy of the predicting results. The model updating based on critical speed of the rotor takes the residual of simulation and test data as the objective value. In this study, the reference model is adopted to represent the test model. By modifying the design parameters in the simplified model using the model updating technique, the 1D model can accurately represent the dynamic characteristics of the reference/test model. The core of the model updating is an iterative optimization problem. The objective function and constraints of the model updating method based on parameter sensitivity are, respectively, defined as follows: where f e and f p ðpÞ are the experimental and analytical where p 0 is the initial values of design parameters. S is the sensitivity matrix of the design parameters. Δp = p − p 0 represents the error of parameters. Thus, the model updating problem can be solved by solving the following equation: This is a set of linear equations that can be solved by an iterative optimization process. In this paper, design variables are screened by sensitivity calculation. The simulated results from the finite element calculation of a 3D model are adopted instead of the test values in model updating. It is noted that when the test values are available, the proposed method can directly be adopted to update the initial simplified model.

Results and Discussion
3.1. Finite Element Modeling. As a reference model for the model updating, a 3D finite element model is established as shown in Figure 5 with 3300 hexahedral elements. The connections in the rotor are considered as integrated. The material parameters for the simulated 3D model are shown in Table 1. The first three free modes are obtained from ANSYS. The mode shapes and natural frequencies are shown in Figures 6(a)-6(c) and Table 2.

Validation of the Three-Dimensional Model.
In order to provide an accurate 3D model in the subsequent model updating process, a modal test as shown in Figure 7 is used to verify the accuracy of the model. The pretightening force of the central rod is 1 × 10 5 N. The elastic rope suspension is used to simulate the free-free boundary. The acceleration sensors are located far away from the joint line of the mode which is shown in Figure 7(c). Four points are hammered circumferentially in a cross-section, and twenty-one points are hammered axially in a generating line. The total number of hammer points is eighty-four. After the data acquisition, the two-dimensional mode shapes are chosen for comparison as shown in Figures 6(d)-6(f). The damping coefficients obtained from the modal test are α = −9:4 × 10 −7 , and β = 7:11. This set of damping coefficients is also applied to the simulation. The mode shapes obtained by the modal test and simulation are shown in Figure 6. The natural frequencies are shown in Table 2. It can be seen from the comparison results that the first three-order bending modes and frequencies obtained from the simulation are consistent with those obtained from the experiment. Therefore, an accurate 3D model of the central rod rotor is obtained, and then, the discrete mass model proposed in this paper will be modified with the critical speed calculated by this model as the target.

Critical Speed Analysis.
The initial values of parameters required for the modeling of the discrete mass model are listed in Table 3 and Table 4. The stiffness of the two elastic supports of the rotor is k a = 1:5 × 10 7 N/m and k b = 1:0 × 10 7 N/m, respectively. Modal analysis on the initial discrete mass model and the 3D model of the rotor is, respectively, conducted when the pretightening force is 100000 N. Because the disk-drum system has the main influence on the vibration of the rotor system, the local vibration of the tie rod is not considered. The mode shapes of two supporting points and five disks are used as the mode shapes of the rotor. The corresponding mode shapes of the first three critical speeds for the two models are shown in Figure 8. The mode shapes corresponding to the first three critical speeds of the rotor obtained by the proposed method are demonstrated in (a), (b), and (c), respectively. (d), (e), and (f) are the first three mode shapes calculated using the 3D finite element model in ANSYS. The first three modes represent the translational, pitch, and bending mode, respectively. As can be seen, the mode shape obtained from the discrete mass model is in good agreement with that obtained from the 3D model. The first three critical speeds and the errors between the two models are listed in Table 5. The errors in critical speed due to the

Sensitivity Analysis and Updating Parameter Selection.
The error of the second and third critical speed of the 1D model is not small. From the results of modal analysis, it can be seen that the modes corresponding to the second and third critical speed are pitching and bending modes, respectively. The pitch mode and bending mode are mainly affected by the support stiffness (k a , k b ) and the bending stiffness (E 1 I 1 , E 2 I 2 , E 3 I 3 , E 4 I 4 , E 5 I 5 , E 6 I 6 ) of the shaft segment, respectively. Therefore, the support stiffness and the bending stiffness of the shaft segment are selected as the parameters to be updated. The sensitivity of critical speed to support stiffness is shown in Figure 9. Results show that the first parameter (left support stiffness) is more sensitive to the second critical speed. The sensitivity of the critical speed to the bending stiffness of six sections of drum shaft is shown in Figure 10. It is obvious that the bending stiffness of the first and second segments has a great influence on the third critical speed. Since the magnitude of the critical speed sensitivity to these two kinds of parameters is not in the same order, they are expressed separately. Based on the above sensitivity analysis, the left support stiffness k a and the flexural stiffness of the first two axle segments denoted as E 1 I 1 and E 2 I 2 are finally selected as the parameters to be updated in modal updating.

Model
Updating. According to the corrected parameters that have been selected, take the minimum error of the first three order critical speed as the optimization target. The convergence of the error in critical speeds verses the number of iterations is shown in Figure 11. The result shows that the    9 International Journal of Aerospace Engineering error of critical speed decreases rapidly at the beginning and converges to a smaller value at the later stage as the number of iterations increases. The variation of parameters with the number of iterations is shown in Figure 12. The errors in critical speed before and after model updating are compared in Table 6. Results show that errors in all the first three critical speeds are reduced to less than 1% after model modification.
3.6. Prediction of Response. In order to verify the validity of the discrete mass model proposed in this paper, the measured damping coefficients are substituted into the updated rod rotor for unbalance response analysis. Limited by the test conditions, the unbalanced response test cannot be carried out temporarily. Therefore, the validated 3D model is used as the reference model for response calculation. The unbalances are 6 × 10 −6 kg·m at the first compressor disk and 8 × 10 −6 kg·m at the gas turbine disk. The tension of the central rod is 1 × 10 5 N. The unbalanced response before and after modal updating and the results obtained by the 3D model are shown in Figure 13. As can be seen, the prediction of the updated model is consistent with that of the 3D model. The accuracy of the proposed method is verified.

Conclusions
A model updating technique based on the discrete mass model with critical speed as the updating target is proposed, which is applied to a central tie rod rotor with preload. The Lagrange equation is used to derive the governing equation of the system. A 3D rotor validated by modal tests is used for critical speed analysis, as a reference for 1D model updating. The correction parameters were selected by the sensitivity method. The accuracy of the updated model is verified by the unbalanced response.
Simulation results show that the accuracy of the model can be effectively improved by using the model updating technique. The error of the first three critical speeds is reduced from more than 5% to less than 1%. The updated model is verified for practicality by unbalanced response using damping coefficient measured by the test. Results show that the method presented in this research work can be used to simulate a complex preloaded rotor system with high efficiency and accuracy.
Although the method in this paper improves the computational efficiency by reducing the degree of freedom, this technique can only be applied to the preliminary qualitative understanding of dynamic characteristics at the rotor design stage. In the future, the model updating based on the test response should be carried out, which has more guiding significance for the identification of supporting parameters and connection parameters.

Data Availability
The data used to support the findings of this study are included within the article.

Conflicts of Interest
The authors declare there is no conflict of interest regarding the publication of this paper.

10
International Journal of Aerospace Engineering