# An Efficient Dynamic Modeling Technique for a Central Tie Rod Rotor

Authors: Jingze Liu; Qingguo Fei; Shaoqing Wu; Zhenhuan Tang; Sanfeng Liao; Dahai Zhang
Journal: International Journal of Aerospace Engineering
 (2021)

---

## Abstract

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 one-dimensional 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 one-dimensional 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.


---

## Body

## 1. 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 Figure1. 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].Figure 1
Cross-section of a typical central tie rod rotor.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 time-consuming. 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 Section2. 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, 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.Figure 2
Flowchart of the methodology.


## 2. Methodology

### 2.1. 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 in Figure3. m1, m2, m3, m4, m5, m6, and m12 represent the lumped mass of the disk drum. m1, m7, m8, m9, m10, m11, and m12 represent the lumped mass of the central tie rod. The disk drum and central tie rod share node 1 and node 12. ka and kb represent the support stiffness of the system.Figure 3
Mass distribution of a central tie rod rotor.


### 2.2. 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 iska and kb, respectively. The distance between the ith disk and the left fulcrum is ai. The distance between the ith mass point on the central tie rod and the left fulcrum is bi. The span between two elastic supports is l. The mass, polar moment of inertia, and equatorial moment of inertia of each point are mi, Jpi, and Jdi, respectively. Each lumped mass point contains four degrees of freedom. xi and yi are the translations along the x- and y-direction at the ith point, respectively. θxi and θyi are the rotations along the x- and y-direction at the ith 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 rotorTr can be expressed as
(1)Tr=12∑i=1nmiẋi2+ẏi2+12∑i=1nJdiθ̇xi2+θ̇yi2+12∑i=1nJpiΩ2−∑i=1nJpiΩθ̇yiθxi.The potential energy of the rotorVr can be expressed as
(2)Vr=12qTKrq,where
(3)q=x1,θy1,⋯,xn,θyn,y1,θx1,⋯,yn,θxnT,(4)Kr=K00K.SymbolK 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 Figure4, assume that the microsegment of the bar is subjected to end forces and external load. According to the shear equilibrium equation, we have
(5)Fsx−Fsx+dFsx+qxdx=0,where Fsx is the shear force on the cross-section. qx is the distribution force per unit length.Figure 4
Force sketch of compressed microsegment.Ifqx equals to zero, then
(6)dFsxdx=qx=0,(7)Mx+dMx−FNwx−Fsx+dFsxdx−Mx=0,where Mx is the moment on the cross-section. FN is the axial force. wx is the deflection of the rod end.From Equation (7), the following equation can be obtained:
(8)dMx−FNwx−Fsx+dFsxdx=0.Eliminating high order small quantity, taking the second derivative ofx, we have
(9)d2Mxdx2−dFsxdx−FNd2wxdx2=0.For straight beams with equal cross-section, the relationship between bending moment and deformation is(10)EIw″=−Mx,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
(11)w4+FNEIw″=0.The general solution of Equation (11) is
(12)w=C1cosFNEIx+C2sinFNEIx+C3x+C4.For the stretched segment, changing theFN in Equation (7) into -FN, the deductions are the same as those for the compressed microsegment. The general solution of the deflection for the stretched segment is as follows:
(13)w=C1eFN/EIx+C2eFN/EIx+C3x+C4.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:(14)ddt∂T∂q̇j−∂T∂qj+∂V∂qj=Qj.By taking the generalized force as zero, the steady-state eddy differential equation of the rotor can be obtained:(15)MMu¨1u¨2+ΩJ−Ju̇1u̇2+KKu1u2=00,where
(16)M=m1Jd1⋱mnJdn,wheremn is the mass of the nth mass point, and Jdn is the diameter moment of inertia of the nth mass point.
(17)J=0Jp1⋱0Jpn,whereJpn is the polar moment of inertia at the nth mass point.The displacement in thexoz plane can be expressed as
(18)u1=x1,θy1,⋯,xn,θynT.The displacement in theyoz plane can be expressed as
(19)u2=y1,θx1,⋯,yn,θxnT.The system damping is assumed to be proportional damping, so the final governing equation can be obtained.(20)Mru¨+Cr−ωGru̇+Kru=Fr,whereMr is the mass matrix, Gr is the gyroscopic matrix, Kr is the stiffness matrix, and Fr is the generalized force vector. The damping of the system Cr is considered as proportional damping, Cr=αKr+βMr, 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.


### 2.3. 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:(21)minWfRp2,Rp=fe−fpps.t.p1≤p≤p2,wherefe and fpp are the experimental and analytical values of critical speed, respectively. Rp is a residual term. p1 and p2 represent the lower and upper bounds of the design parameters, respectively. Wf represents the weighted matrix for all the eigenvalues.Performing the first-order Taylor expansion onfpp at the initial design point, we have
(22)fpp=fpp0+SΔp,where p0 is the initial values of design parameters. S is the sensitivity matrix of the design parameters. ∆p=p−p0 represents the error of parameters. Thus, the model updating problem can be solved by solving the following equation:
(23)WfSΔp=Wffe−fpp0.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.


## 2.1. 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 in Figure3. m1, m2, m3, m4, m5, m6, and m12 represent the lumped mass of the disk drum. m1, m7, m8, m9, m10, m11, and m12 represent the lumped mass of the central tie rod. The disk drum and central tie rod share node 1 and node 12. ka and kb represent the support stiffness of the system.Figure 3
Mass distribution of a central tie rod rotor.


## 2.2. 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 iska and kb, respectively. The distance between the ith disk and the left fulcrum is ai. The distance between the ith mass point on the central tie rod and the left fulcrum is bi. The span between two elastic supports is l. The mass, polar moment of inertia, and equatorial moment of inertia of each point are mi, Jpi, and Jdi, respectively. Each lumped mass point contains four degrees of freedom. xi and yi are the translations along the x- and y-direction at the ith point, respectively. θxi and θyi are the rotations along the x- and y-direction at the ith 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 rotorTr can be expressed as
(1)Tr=12∑i=1nmiẋi2+ẏi2+12∑i=1nJdiθ̇xi2+θ̇yi2+12∑i=1nJpiΩ2−∑i=1nJpiΩθ̇yiθxi.The potential energy of the rotorVr can be expressed as
(2)Vr=12qTKrq,where
(3)q=x1,θy1,⋯,xn,θyn,y1,θx1,⋯,yn,θxnT,(4)Kr=K00K.SymbolK 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 Figure4, assume that the microsegment of the bar is subjected to end forces and external load. According to the shear equilibrium equation, we have
(5)Fsx−Fsx+dFsx+qxdx=0,where Fsx is the shear force on the cross-section. qx is the distribution force per unit length.Figure 4
Force sketch of compressed microsegment.Ifqx equals to zero, then
(6)dFsxdx=qx=0,(7)Mx+dMx−FNwx−Fsx+dFsxdx−Mx=0,where Mx is the moment on the cross-section. FN is the axial force. wx is the deflection of the rod end.From Equation (7), the following equation can be obtained:
(8)dMx−FNwx−Fsx+dFsxdx=0.Eliminating high order small quantity, taking the second derivative ofx, we have
(9)d2Mxdx2−dFsxdx−FNd2wxdx2=0.For straight beams with equal cross-section, the relationship between bending moment and deformation is(10)EIw″=−Mx,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
(11)w4+FNEIw″=0.The general solution of Equation (11) is
(12)w=C1cosFNEIx+C2sinFNEIx+C3x+C4.For the stretched segment, changing theFN in Equation (7) into -FN, the deductions are the same as those for the compressed microsegment. The general solution of the deflection for the stretched segment is as follows:
(13)w=C1eFN/EIx+C2eFN/EIx+C3x+C4.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:(14)ddt∂T∂q̇j−∂T∂qj+∂V∂qj=Qj.By taking the generalized force as zero, the steady-state eddy differential equation of the rotor can be obtained:(15)MMu¨1u¨2+ΩJ−Ju̇1u̇2+KKu1u2=00,where
(16)M=m1Jd1⋱mnJdn,wheremn is the mass of the nth mass point, and Jdn is the diameter moment of inertia of the nth mass point.
(17)J=0Jp1⋱0Jpn,whereJpn is the polar moment of inertia at the nth mass point.The displacement in thexoz plane can be expressed as
(18)u1=x1,θy1,⋯,xn,θynT.The displacement in theyoz plane can be expressed as
(19)u2=y1,θx1,⋯,yn,θxnT.The system damping is assumed to be proportional damping, so the final governing equation can be obtained.(20)Mru¨+Cr−ωGru̇+Kru=Fr,whereMr is the mass matrix, Gr is the gyroscopic matrix, Kr is the stiffness matrix, and Fr is the generalized force vector. The damping of the system Cr is considered as proportional damping, Cr=αKr+βMr, 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.


## 2.3. 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:(21)minWfRp2,Rp=fe−fpps.t.p1≤p≤p2,wherefe and fpp are the experimental and analytical values of critical speed, respectively. Rp is a residual term. p1 and p2 represent the lower and upper bounds of the design parameters, respectively. Wf represents the weighted matrix for all the eigenvalues.Performing the first-order Taylor expansion onfpp at the initial design point, we have
(22)fpp=fpp0+SΔp,where p0 is the initial values of design parameters. S is the sensitivity matrix of the design parameters. ∆p=p−p0 represents the error of parameters. Thus, the model updating problem can be solved by solving the following equation:
(23)WfSΔp=Wffe−fpp0.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.


## 3. 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 Figure5 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.Figure 5
3D finite element model.Table 1
Material parameters for the 3D model.
Elastic modulus (GPa)Density (kg/m3)First compressor disk1967800Second compressor disk1214480Third compressor disk1214480Centrifugal impeller disk1214480Gas turbine disk2148200Central tie rod2047850Figure 6
Experimental and simulated modes of vibration: (a), (b), and (c) are from simulated; (d), (e), and (f) are from experiment.
(a)
Mode 1(b)
Mode 2(c)
Mode 3(d)
Mode 1(e)
Mode 2(f)
Mode 3Table 2
Errors of frequencies (Hz).
SimulatedExperimentalErrors (%)Mode 1127512551.6Mode 2255625251.2Mode 3364835841.8


### 3.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 Figure7 is used to verify the accuracy of the model. The pretightening force of the central rod is 1×105N. 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.Figure 7
Modal test of the rod rotor: (a) rod rotor; (b) data acquisition system; (c) schematic diagram of experimental.
(a)(b)(c)


### 3.3. Critical Speed Analysis

The initial values of parameters required for the modeling of the discrete mass model are listed in Table3 and Table 4. The stiffness of the two elastic supports of the rotor is ka=1.5×107N/m and kb=1.0×107N/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 simplification process can be reduced via the model updating.Table 3
Modeling parameters of shaft of the 1D model.
ParameterValue between theith node and the jth nodeLength (mm)L12L23L34L45L56L6127575559010565L17L78L89L910L1011L1112658585858560Outer diameter (mm)D12D23D34D45D56D61250751101109050D17D78D89D910D1011D1112404040404040Inner diameter (mm)d12d23d34d45d56d61240651001008040d17d78d89d910d1011d1112323232323232Elastic modulus (GPa)E12E23E34E45E56E612196121121121214214E17E78E89E910E1011E1112204204204204204204Table 4
Modeling parameters of lumped mass for the 1D model.
Mass (kg)Diameter moment of inertia (10-4 kg·m2)Polar moment of inertia (10-4 kg·m2)Distance of centroid (mm)m10.4Jd11.6Jp11.8a171.5m21.2Jd210.6Jp214.2a2149.1m30.8Jd39.0Jp315.1a3209.3m40.9Jd412.5Jp420.4a4296.2m55.5Jd5199.9Jp5379.3a5398.4m65.1Jd698.8Jp6178.7a6465m70.3Jd72.4Jp71.0b165m80.3Jd82.4Jp81.0b2150m90.3Jd92.4Jp91.0b3235m100.3Jd102.4Jp101.0b4320m110.3Jd112.4Jp111.0b5405m120.3Jd121.1Jp121.4Figure 8
Mode shapes of the 1D model and 3D model: (a), (b), and (c) are from the 1D model; (d), (e), and (f) are from the 3D model.
(a)
Mode 1(b)
Mode 2(c)
Mode 3(d)
Mode 1(e)
Mode 2(f)
Mode 3Table 5
Errors of critical speed.
Mode1D (rpm)3D model (rpm)Errors (%)11141611460-0.38227105255006.293109668999659.71


### 3.4. 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 (ka, kb) and the bending stiffness (E1I1, E2I2, E3I3, E4I4, E5I5, E6I6) 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 Figure9. 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 ka and the flexural stiffness of the first two axle segments denoted as E1I1 and E2I2 are finally selected as the parameters to be updated in modal updating.Figure 9
Sensitivity of critical speed to supporting stiffness.Figure 10
Sensitivity of critical speed to bending stiffness.


### 3.5. 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 Figure11. The result shows that the 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.Figure 11
Convergence of the percentage errors in critical speeds verses the number of iterations.Figure 12
Variation of parameters with the number of iterations.Table 6
Errors in critical speeds after model updating.
Mode1D (rpm)3D model (rpm)Errors (%)11142711460-0.29225540255000.163100514999650.55


### 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 are6×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×105N. 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.Figure 13
Comparisons of unbalanced response.
(a)
First compressor disk(b)
Gas turbine disk


## 3.1. Finite Element Modeling

As a reference model for the model updating, a 3D finite element model is established as shown in Figure5 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.Figure 5
3D finite element model.Table 1
Material parameters for the 3D model.
Elastic modulus (GPa)Density (kg/m3)First compressor disk1967800Second compressor disk1214480Third compressor disk1214480Centrifugal impeller disk1214480Gas turbine disk2148200Central tie rod2047850Figure 6
Experimental and simulated modes of vibration: (a), (b), and (c) are from simulated; (d), (e), and (f) are from experiment.
(a)
Mode 1(b)
Mode 2(c)
Mode 3(d)
Mode 1(e)
Mode 2(f)
Mode 3Table 2
Errors of frequencies (Hz).
SimulatedExperimentalErrors (%)Mode 1127512551.6Mode 2255625251.2Mode 3364835841.8


## 3.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 Figure7 is used to verify the accuracy of the model. The pretightening force of the central rod is 1×105N. 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.Figure 7
Modal test of the rod rotor: (a) rod rotor; (b) data acquisition system; (c) schematic diagram of experimental.
(a)(b)(c)


## 3.3. Critical Speed Analysis

The initial values of parameters required for the modeling of the discrete mass model are listed in Table3 and Table 4. The stiffness of the two elastic supports of the rotor is ka=1.5×107N/m and kb=1.0×107N/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 simplification process can be reduced via the model updating.Table 3
Modeling parameters of shaft of the 1D model.
ParameterValue between theith node and the jth nodeLength (mm)L12L23L34L45L56L6127575559010565L17L78L89L910L1011L1112658585858560Outer diameter (mm)D12D23D34D45D56D61250751101109050D17D78D89D910D1011D1112404040404040Inner diameter (mm)d12d23d34d45d56d61240651001008040d17d78d89d910d1011d1112323232323232Elastic modulus (GPa)E12E23E34E45E56E612196121121121214214E17E78E89E910E1011E1112204204204204204204Table 4
Modeling parameters of lumped mass for the 1D model.
Mass (kg)Diameter moment of inertia (10-4 kg·m2)Polar moment of inertia (10-4 kg·m2)Distance of centroid (mm)m10.4Jd11.6Jp11.8a171.5m21.2Jd210.6Jp214.2a2149.1m30.8Jd39.0Jp315.1a3209.3m40.9Jd412.5Jp420.4a4296.2m55.5Jd5199.9Jp5379.3a5398.4m65.1Jd698.8Jp6178.7a6465m70.3Jd72.4Jp71.0b165m80.3Jd82.4Jp81.0b2150m90.3Jd92.4Jp91.0b3235m100.3Jd102.4Jp101.0b4320m110.3Jd112.4Jp111.0b5405m120.3Jd121.1Jp121.4Figure 8
Mode shapes of the 1D model and 3D model: (a), (b), and (c) are from the 1D model; (d), (e), and (f) are from the 3D model.
(a)
Mode 1(b)
Mode 2(c)
Mode 3(d)
Mode 1(e)
Mode 2(f)
Mode 3Table 5
Errors of critical speed.
Mode1D (rpm)3D model (rpm)Errors (%)11141611460-0.38227105255006.293109668999659.71


## 3.4. 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 (ka, kb) and the bending stiffness (E1I1, E2I2, E3I3, E4I4, E5I5, E6I6) 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 Figure9. 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 ka and the flexural stiffness of the first two axle segments denoted as E1I1 and E2I2 are finally selected as the parameters to be updated in modal updating.Figure 9
Sensitivity of critical speed to supporting stiffness.Figure 10
Sensitivity of critical speed to bending stiffness.


## 3.5. 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 Figure11. The result shows that the 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.Figure 11
Convergence of the percentage errors in critical speeds verses the number of iterations.Figure 12
Variation of parameters with the number of iterations.Table 6
Errors in critical speeds after model updating.
Mode1D (rpm)3D model (rpm)Errors (%)11142711460-0.29225540255000.163100514999650.55


## 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 are6×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×105N. 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.Figure 13
Comparisons of unbalanced response.
(a)
First compressor disk(b)
Gas turbine disk


## 4. 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.



---

*Source: 6618828-2021-06-08.xml*