# Dynamic Analysis of Flexible Rotor Based on Transfer Symplectic Matrix

Authors: Hao Deng; Xi Fang; Huachun Wu; Yiming Ding; Jinghu Yu; Xiaoyu Zhang; Hong Li; Chenwei Liu
Journal: Shock and Vibration
 (2019)

---

## Abstract

In view of the numerical instability and low accuracy of the traditional transfer matrix method in solving the high-order critical speed of the rotor system, a new idea of incorporating the finite element method into the transfer matrix is proposed. Based on the variational principle, the transfer symplectic matrix of gyro rotors suitable for all kinds of boundary conditions and supporting conditions under the Hamilton system is derived by introducing dual variables. To verify the proposed method in rotor critical speed, a numerical analysis is adopted. The simulation experiment results show that, in the calculation of high-order critical speed, especially when exceeding the sixth critical speed, the numerical accuracy of the transfer symplectic matrix method is obviously better than that of the reference method. The relative errors between the numerical solution and the exact solution are 0.0347% and 0.2228%, respectively, at the sixth critical speed. The numerical example indicates the feasibility and superiority of the method, which provides the basis for the optimal design of the rotor system.


---

## Body

## 1. Introduction

With the development of modern industrial technology, the application prospect and running speed of rotary machinery have been greatly improved. Correspondingly, the performance and design requirements of rotary machinery have also been put forward with greater challenges [1]. In order to make the rotor work smoothly, avoid resonance and other vibration faults, and ensure high rotation accuracy, the study of the natural frequency and mode of the rotor dynamics system, the critical speed, and the unbalanced response has become an indispensable important content of rotor dynamics. Critical speed and vibration mode, which are important dynamic characteristic parameters of rotating machinery, are closely related to various structural parameters of the rotor system itself, such as the influence of unbalance and support stiffness on critical speed. Therefore, it is of great significance to analyze the sensitivity of the dynamic characteristic parameters of the rotor system to the structural parameters of the system and to optimize the rotor structure and shorten the design period [2].So far, the numerical analysis methods of rotor dynamics are mainly divided into two categories: transfer matrix method and finite element method. The finite element method combines the variational principle of mathematics with the energy principle of mechanics [3, 4]. According to the principle of minimum potential energy of elasticity, when a deformed elastomer is in an equilibrium state under the action of external forces, the curve that minimizes the total potential energy is the real deformation curve among many possible deformation permissible curves. The expression of the rotor finite element method is concise and standard, which has prominent advantages when solving the problems of the large complex system composed of rotor and surrounding structures, but this method will increase the system freedom degree and result in a lot of time and space spending [5–7].For the multisupported multidisk rotor system, when the system has more degrees of freedom, whether by establishing the motion differential equation to solve the critical speed and the unbalanced response of the system or by solving the eigenvalue problem, the calculation workload is very large. Therefore, a new computational method, the classical transfer matrix method (Prohl transfer matrix method), was developed. This method is widely used because the order of the matrix will not increase with the increase of the degrees of freedom of the system and less storage units are needed [8–10]. Nevertheless, with the increase of frequency of trial computation, the accuracy decreases, especially in solving the critical speeds and unbalance response of the rotor structure system, which is of large scale and high speed. The vibration mode of the structural system firstly appears numerical instability, which is mainly manifested as a sharp increase in the amplitude at the end of the shaft system, and then the root of the frequency equation begins to be errors and omissions, resulting in a larger error in the calculation results. The Riccati transfer matrix method is then proposed to transform the original two-point boundary value problem of the differential equation into a one-point initial value problem through the Riccati transformation, which improves the numerical stability and calculation accuracy while retaining the advantages of the classical transfer matrix method [11–14]. However, when solving the critical speed, there may be some problems such as increasing roots or missing roots because the interval of residual quantity with different signs does not correspond to the interval of residual quantity root.The transfer matrix method mentioned above is to analyze the dynamic characteristics of the rotor under the Lagrange system. Academician Zhong Wanxie introduced some important issues of classical mechanics into the Hamilton system; he studied the dynamics of the gyro rotor from a new perspective and proved the effectiveness of analyzing the gyro rotor system under the Hamilton system [15–17].Consider the advantages and disadvantages of each method mentioned above, in order to make the rotor work smoothly, avoid resonance, and guarantee its high rotation accuracy; the idea of finite element is incorporated into the transfer matrix in this paper. Starting from the variational principle, dual variables are introduced to deduce the transfer matrix of the gyro rotor under the Hamilton system. Then, numerical examples are given to verify the feasibility and superiority of this method. Since the transfer matrix between state vectors is derived by symplectic method, it has significant advantages in long-term numerical stability. In addition, compared with the traditional transfer matrix method, this method avoids a series of complex operations such as solving the aggregation model of the beam, which is easy to understand, and can obtain relatively accurate results when calculating the high-order vortex frequency.A brief outline of this paper is as follows. In the second section, we firstly review the classical transfer matrix methods such as Prohl transfer matrix method and Riccati transfer matrix method and then introduce the dual-rotor system and its transfer matrix method. In the third section, we derive the transfer symplectic matrix between state vectors of the rotor system under the Hamiltonian system from the variational principle and briefly explained that this type of transfer matrix is symplectic. In the fourth section, the numerical example is used for verification and different numerical algorithms are compared. In the fifth section, we mainly introduce the conclusions and innovations.


## 2. Theoretical Basis

### 2.1. Classical Transfer Matrix Method

As one of the analysis methods of rotor dynamics, the transfer matrix method is usually used to solve rotor dynamics problems in the frequency domain. As the transfer matrix method does not require the storage and the dimension does not increase with the number of elements (for it is replaced with more matrix multiplications), this kind of method is especially suitable for chain systems such as rotor systems.Figure1 depicts the rotor structure adopted in this section. Before analyzing the dynamics of the subsystem by the transfer matrix method, the rotor subsystem is essential to be discretized and simplified into the lumped mass model [18] (shown in Figure 2). The assembled parts on the rotor, such as the front and rear lamination stacks, are modeled as rigid disks, with consideration of the gyroscopic moment. Considering the x−O−z plane of the discrete rotor model, a typical unit consists of a massless shaft segment and a lumped mass disc [19].Figure 1
Rotor model and structure drawing.
(a)
(b)Figure 2
Shaft torsional vibration system.Considering thei−th element, defining Z as the stator vector containing deflection y, with θ as the slope, M as the bending moment, and Q as the shear force, we obtain(1)Zi=y,θ,M,QiT.For any axial segmenti, there is a certain relationship between the state vectors ZiandZi+1, that is(2)Zi+1=TiZi,where Ti is called the transfer matrix of components between the two sections and the matrix elements can be obtained by force analysis and deformation relation of components.In order to facilitate the study of the relationship between adjacent nodes and axis segments, this paper did the simplification as follows: The bending moment and shear force to the right side of nodei−1 represent the bending moment and shear force to the left side of axis segment i, respectively, while the bending moment and shear force to the left of node i represent the bending moment and shear force to the right side of axis segment i, respectively.According to the bending deformation formula and displacement analysis of beam, the following equation is obtained [20–22]:(3)Zi+1L=AiZiR,where the field transfer matrix Ai is(4)Ai=1lili22EiIili31−υ6EiIi01liEiIili22EiIi001li0001,in which υ=ks6EiIi/GiAili2 represents the shear effect coefficient, ks is called the cross-sectional shape coefficient, and ksli/GiAi is the influence of shear deformation on bending deformation.For isotropic rotor systems, when the system moves in a simple harmonic motion with the angular velocityω, we get y¨i=−ω2yi. Suppose that the j−th node has a spring bearing and the spring stiffness is kj. According to rigid body dynamics and D’Alembert’s principle, the relationship between the state vectors at the left and right ends of the rotating disk supported by elasticity can be obtained as [23](5)ZiR=CiZiL,where the transfer matrix Ci is(6)Ci=100001000−Jdi−JpiΩωω210miω2−kj001,in which Jdi is the moment of inertia of diameter, Jpi is the polar moment of inertia, Ω is the angular velocity of axis rotation, and ω is the angular speed of revolution.By combining the flexible rotor and rigid disk with a transfer element, the relation of the state vectors fromi−th to i+1−th can be written as(7)Zi+1L=AiZiR=AiCiZiL≜TiZiL.When the transfer matrix of each part of the rotor system is obtained from the above calculation, the whole system is constructed by successively multiplying the matrix of each part according to equation (7):(8)ZN+1=TNZN=⋯=TNTN−1⋯T2T1Z1=DZ1,where the assembled transferring matrix D is(9)D=d11d12d13d14d21d22d23d24d31d32d33d34d41d42d43d44.It is found that the value of the element containingω in the total transfer matrix D is only restricted by the boundary conditions on both sides, so the critical angular velocity of the system can be obtained by solving the ω, which satisfies the boundary conditions [24].When the traditional transfer matrix method is used to calculate the response of high-frequency dynamics or systems with an excessive number of nodes, the numerical instability occurs. In order to solve the problem, the Riccati transfer matrix method appears. The Riccati transfer matrix method occurs which converts a numerically unstable two-point boundary value problem into a numerically stable value problem. Through the Riccati transform, the state vectorZi+1 is converted into(10)Zi+1=fei+1=u11u12u21u22fei,where f and e represent the force part of the state vector and the displacement part of the state vector, respectively. A generalized Riccati transformation at the i‐th section is given as follows:(11)fi=siei,where si is the so-called Riccati transfer matrix. And through the derivation, we get the recursion formula(12)si+1=u11s+u12iu21s+u22i−1.


### 2.2. Dual-Rotor System Transfer Matrix Method

Consider a dual-rotor system with two intermediate support, as shown in Figure3. Defining the state vectors of the i‐th element of the two rotors as Zi1andZi2, respectively:(13)Zi1=y1,θ1,M1,Q1iT,Zi2=y2,θ2,M2,Q2iT.Figure 3
Schematic diagram of dual-rotor system.As for the noncoupling element, according to the transfer matrix method of single rotor system, their respective transfer relations can be obtained:(14)Z1Z2i+1=Ti100Ti2Z1Z2i,denoted by(15)Zi+1=TiZi,Ti1=1+l36EIC1C2l+l22EIC3l22EIl36EIC1l22EIC21+lEIC3lEIl22EIlC2C31lC2001i1,Ti2=1+l36EIC1C2l+l22EIC3l22EIl36EIC1l22EIC21+lEIC3lEIl22EIlC2C31lC2001i2,in which Ti1andTi2 are the transfer matrix of the two rotors. Their form can be derived from the Prohl, Riccati, or other transfer matrix methods. C1=1−v, C2=mω2−k, and C3=−Jd−JpΩ/ωω2+K, where kandK are defined as 0 when there is no spring bearing in the i‐th element.As for the coupling element, assume that the coupling point is in thei‐th unit. According to the balance of forces, the coordination equations are obtained:(16)y1θ1M1Q1iR=G1G2y1θ1M1Q1y2θ2M2Q2iLT,y2θ2M2Q2iR=G2G1y1θ1M1Q1y2θ2M2Q2iLT,so that we can get state vectors on both sides of the coupling element:(17)ZiR=TGiZiL,where TGi represents the transfer matrix for the coupling element, denoted as(18)TGi=G1G2G2G1i,G1=11Kh1−k1i,G2=−Khki,where k stands for support stiffness of the spring and Kh is the torque stiffness.


### 2.3. Critical Speed and Mode of Vibration

Supposing that the assembled transferring matrix is denoted byZN+1R=DZ1R, according to the boundary conditions of the model(19)M1=0,Q1=0,MN+1=0,QN+1=0.Combine equations (17) and (19), there is(20)d31d32d42d42yθ0=00.Equation (20) is a system of homogeneous algebraic equations, where each element is a function of vortex frequency, and its coefficient determinant Δω2 called surplus is(21)Δω2=d31d32d42d42.The condition for the existence of nonzero solutions to the homogeneous algebraic equations is that the determinant of coefficients is equal to zero, thus obtaining the rotor frequency equation:(22)Δω2=Δd31d32d42d42=0.The valueω satisfying equation (22) is not only a synchronous forward vortex angular velocity of the rotor but also the critical angular velocity of the rotor.


## 2.1. Classical Transfer Matrix Method

As one of the analysis methods of rotor dynamics, the transfer matrix method is usually used to solve rotor dynamics problems in the frequency domain. As the transfer matrix method does not require the storage and the dimension does not increase with the number of elements (for it is replaced with more matrix multiplications), this kind of method is especially suitable for chain systems such as rotor systems.Figure1 depicts the rotor structure adopted in this section. Before analyzing the dynamics of the subsystem by the transfer matrix method, the rotor subsystem is essential to be discretized and simplified into the lumped mass model [18] (shown in Figure 2). The assembled parts on the rotor, such as the front and rear lamination stacks, are modeled as rigid disks, with consideration of the gyroscopic moment. Considering the x−O−z plane of the discrete rotor model, a typical unit consists of a massless shaft segment and a lumped mass disc [19].Figure 1
Rotor model and structure drawing.
(a)
(b)Figure 2
Shaft torsional vibration system.Considering thei−th element, defining Z as the stator vector containing deflection y, with θ as the slope, M as the bending moment, and Q as the shear force, we obtain(1)Zi=y,θ,M,QiT.For any axial segmenti, there is a certain relationship between the state vectors ZiandZi+1, that is(2)Zi+1=TiZi,where Ti is called the transfer matrix of components between the two sections and the matrix elements can be obtained by force analysis and deformation relation of components.In order to facilitate the study of the relationship between adjacent nodes and axis segments, this paper did the simplification as follows: The bending moment and shear force to the right side of nodei−1 represent the bending moment and shear force to the left side of axis segment i, respectively, while the bending moment and shear force to the left of node i represent the bending moment and shear force to the right side of axis segment i, respectively.According to the bending deformation formula and displacement analysis of beam, the following equation is obtained [20–22]:(3)Zi+1L=AiZiR,where the field transfer matrix Ai is(4)Ai=1lili22EiIili31−υ6EiIi01liEiIili22EiIi001li0001,in which υ=ks6EiIi/GiAili2 represents the shear effect coefficient, ks is called the cross-sectional shape coefficient, and ksli/GiAi is the influence of shear deformation on bending deformation.For isotropic rotor systems, when the system moves in a simple harmonic motion with the angular velocityω, we get y¨i=−ω2yi. Suppose that the j−th node has a spring bearing and the spring stiffness is kj. According to rigid body dynamics and D’Alembert’s principle, the relationship between the state vectors at the left and right ends of the rotating disk supported by elasticity can be obtained as [23](5)ZiR=CiZiL,where the transfer matrix Ci is(6)Ci=100001000−Jdi−JpiΩωω210miω2−kj001,in which Jdi is the moment of inertia of diameter, Jpi is the polar moment of inertia, Ω is the angular velocity of axis rotation, and ω is the angular speed of revolution.By combining the flexible rotor and rigid disk with a transfer element, the relation of the state vectors fromi−th to i+1−th can be written as(7)Zi+1L=AiZiR=AiCiZiL≜TiZiL.When the transfer matrix of each part of the rotor system is obtained from the above calculation, the whole system is constructed by successively multiplying the matrix of each part according to equation (7):(8)ZN+1=TNZN=⋯=TNTN−1⋯T2T1Z1=DZ1,where the assembled transferring matrix D is(9)D=d11d12d13d14d21d22d23d24d31d32d33d34d41d42d43d44.It is found that the value of the element containingω in the total transfer matrix D is only restricted by the boundary conditions on both sides, so the critical angular velocity of the system can be obtained by solving the ω, which satisfies the boundary conditions [24].When the traditional transfer matrix method is used to calculate the response of high-frequency dynamics or systems with an excessive number of nodes, the numerical instability occurs. In order to solve the problem, the Riccati transfer matrix method appears. The Riccati transfer matrix method occurs which converts a numerically unstable two-point boundary value problem into a numerically stable value problem. Through the Riccati transform, the state vectorZi+1 is converted into(10)Zi+1=fei+1=u11u12u21u22fei,where f and e represent the force part of the state vector and the displacement part of the state vector, respectively. A generalized Riccati transformation at the i‐th section is given as follows:(11)fi=siei,where si is the so-called Riccati transfer matrix. And through the derivation, we get the recursion formula(12)si+1=u11s+u12iu21s+u22i−1.


## 2.2. Dual-Rotor System Transfer Matrix Method

Consider a dual-rotor system with two intermediate support, as shown in Figure3. Defining the state vectors of the i‐th element of the two rotors as Zi1andZi2, respectively:(13)Zi1=y1,θ1,M1,Q1iT,Zi2=y2,θ2,M2,Q2iT.Figure 3
Schematic diagram of dual-rotor system.As for the noncoupling element, according to the transfer matrix method of single rotor system, their respective transfer relations can be obtained:(14)Z1Z2i+1=Ti100Ti2Z1Z2i,denoted by(15)Zi+1=TiZi,Ti1=1+l36EIC1C2l+l22EIC3l22EIl36EIC1l22EIC21+lEIC3lEIl22EIlC2C31lC2001i1,Ti2=1+l36EIC1C2l+l22EIC3l22EIl36EIC1l22EIC21+lEIC3lEIl22EIlC2C31lC2001i2,in which Ti1andTi2 are the transfer matrix of the two rotors. Their form can be derived from the Prohl, Riccati, or other transfer matrix methods. C1=1−v, C2=mω2−k, and C3=−Jd−JpΩ/ωω2+K, where kandK are defined as 0 when there is no spring bearing in the i‐th element.As for the coupling element, assume that the coupling point is in thei‐th unit. According to the balance of forces, the coordination equations are obtained:(16)y1θ1M1Q1iR=G1G2y1θ1M1Q1y2θ2M2Q2iLT,y2θ2M2Q2iR=G2G1y1θ1M1Q1y2θ2M2Q2iLT,so that we can get state vectors on both sides of the coupling element:(17)ZiR=TGiZiL,where TGi represents the transfer matrix for the coupling element, denoted as(18)TGi=G1G2G2G1i,G1=11Kh1−k1i,G2=−Khki,where k stands for support stiffness of the spring and Kh is the torque stiffness.


## 2.3. Critical Speed and Mode of Vibration

Supposing that the assembled transferring matrix is denoted byZN+1R=DZ1R, according to the boundary conditions of the model(19)M1=0,Q1=0,MN+1=0,QN+1=0.Combine equations (17) and (19), there is(20)d31d32d42d42yθ0=00.Equation (20) is a system of homogeneous algebraic equations, where each element is a function of vortex frequency, and its coefficient determinant Δω2 called surplus is(21)Δω2=d31d32d42d42.The condition for the existence of nonzero solutions to the homogeneous algebraic equations is that the determinant of coefficients is equal to zero, thus obtaining the rotor frequency equation:(22)Δω2=Δd31d32d42d42=0.The valueω satisfying equation (22) is not only a synchronous forward vortex angular velocity of the rotor but also the critical angular velocity of the rotor.


## 3. Construction of the Transfer Matrix in Symplectic Space

Different from the traditional transfer matrix construction, this paper deduces the transfer symplectic matrix of the linear rotor system in symplectic space from the variational principle. The idea of finite element is incorporated into the transfer matrix, and the transfer symplectic matrix of gyro rotor under the Hamilton system is derived by introducing dual variables from the variational principle. Using this method to calculate the critical speed has the advantages of good numerical stability, high calculation accuracy and easy to compile general calculation program, and it has been verified by numerical examples.Consider the following types of rotor systems: isotropic Timoshenko beams are used at the connection beams. Take a section of beam, which considers the moment of inertia and shear deformation for analysis. Let thez-coordinate be along the axis of the beam before deformation. Let μ be the cross-section shear coefficient, G be the shear modulus, E be the elastic modulus, I be the moment of inertia, A be the cross-section area of the beam, and ux,uy and ψx,ψy be the linear and angular displacements of the cross section, respectively. The deformation energy of this beam can be expressed as(23)V=12∫0lEIdψxdz2+dψydz2+μAGduxdz−ψx2+duydz−ψy2dz.Take thex−O−z perpendicular plane, for example, the internal force of the beam can be expressed as(24)M=EIdψxdz,Q=μAGduxdz−ψx.The dynamic equations of the beam can be deduced to(25)∂Q∂z=ρA∂2ux∂t2,Q−∂M∂z=ρI∂2ψx∂t2.For the sake of analysis, the time domain is converted to the frequency domain and the variable is still represented by the previous notation(26)u=uz,ω⋅e−iωt,ψ=ψz,ω⋅e−iωt,where ω is the rotation angular velocity of the rotor, and then, the dynamic equation (25) can be converted into(27)ddzμAGduxdz−ψx+ρAω2ux=0,ddzEIdψxdz+μAGduxdz−ψx+ρIω2ψx=0.By the same logic, we can get that of they−O−z plane. Thus, we can obtain the kinetic energy of this beam:(28)Ek=12∫0lρAduxdz2+duydz2+ρIdψxdz2+dψydz2+2ΩρIdψxdzψy−dψydzψxdz.Therefore, the Lagrangian density function is(29)L=12ρAduxdz2+duydz2+ρIdψxdz2+dψydz2+2ΩρIdψxdzψy−dψydzψx−EIdψxdz2+dψydz2−μAGduxdz−ψx2+duydz−ψy2.For the convenience of analysis and calculation, the following notation is introduced:q=u,ψT;q˙=u˙,ψ˙T. Take the matrix form(30)L=12qTρAω2ρIω2q−qTΩρIωq−12q˙TμAGEIq˙−q˙T−μAGq−12qTμAGq.Then, the potential energy variation principle can be given for the above equation as follows:(31)δ⋅12∫0lLq,q˙dz=0,Lq,q˙=12q˙TK22q˙+q˙TK21q+12qTK11q.The corresponding variational principle isSqz=∫0LLq,q˙dz,δS=0 and of which(32)K22=μAGEI,K21=−μAG,K11=−ρAω2ρIω2−μAG−2ΩρIω,K12=K21T.Then, the dynamic equation can be re-expressed as(33)K22q¨+K21−K12q˙−K11q=0.Considering that both the variational expression and the dynamic equation only include the displacement type of variables, the method of introducing dual variables is adopted, letp=∂L/∂q˙ and then it can be re-expressed as(34)p=K22q˙+K21q.Substitute equation (34) into the dynamic equation, we can get(35)q˙=Aq+Cp,p˙=Bq−ATp,where A=−K22−1K21,B=K11−K12K22−1K21,andEC=K22−1.Then, the HamiltonianHp,q is introduced as(36)Hp,q=pTAq−12qTBq+12pTCp.The corresponding variational principle isSqz,pz=∫z0zfpTq˙−Hq,pdz,δS=0, and the dual vectors qandp are unified to form the full state vector vz=qT,pTT.The equation above is for the generaln-dimensional displacement vector; in this paper, there is n=2, then p=K22q˙+K21q=μAGu˙−ψ,EIψ˙T=Q.−M. Moreover, at the same time, we can get the Hamiltonian equation of state as(37)v˙=Hv,H=ACB−AT,in which(38)A=11,B=−ρω2A−ρω2I,C=μAG−1EI−1.It is observed that the mechanical meaning of dual vectorp is composed of generalized internal force, namely, shear force Q and bending moment M, which shows the rationality and validity of the introduced dual variables.When solving the ordinary differential equation (37), we let li=Zi+1−Zi and(39)vi+1=eHilivi.Along the direction of theZ-axis, divide the Z-axis into arbitrarily small intervals. The transformation matrix when transferring from section i to section i+1 is(40)Fi=eHili.Now that theli=Zi+1−Zi is small enough, we can expand eHili by Taylor series(41)Fi=eHili=I+liHi+li22!Hi2+li33!Hi3+⋯.Equation (41) indicates that the series converges so that the method put forward is available. Then, vi+1=Fivi, in which Fi is the transfer matrix Ai as mentioned in Section 2.The symplectic algorithm maintains strictly the symplectic structure of the Hamilton system and has unique advantages in long-term numerical stability [25]. Therefore, the transfer symplectic matrix method proposed in this paper has better numerical stability and accuracy in solving the critical speed of the rotor.


## 4. Numerical Results

Traditional single stator and rotor motors, whether dc, asynchronous, or synchronous, have only one mechanical port. In recent years, the concept of dual-rotor motor has been put forward, which enables the design of the motor to realize the dual mechanical shaft and realize the independent transfer of energy on the two mechanical shafts. The dual-rotor motor can greatly reduce the actual size and weight of the motor and greatly improve the working efficiency of the motor, which can well meet the specific requirements of energy saving and speed regulation of the motor, showing better motor performance. Therefore, dual-rotor motor has been widely used in many fields and this paper will take the dual-rotor system as an example for verification analysis.In this section, in order to validate the accuracy and efficiency of the method presented in this paper, a numerical example is given and the critical speeds of each order are compared with those presented in [26, 27] and other numerical algorithms; the numerical analysis and symbolic computations are performed in MATLAB.As an example, we choose the two-rotor system with two intermediate supports, which is shown in Figure4. Some material and geometrical properties are given as E = 206 GPa, μ = 0.3, and ρ = 7850 kg/m3. The detailed original data of the calculation example are shown in Table 1.Figure 4
Computational model of a dual-rotor system.Table 1
Original data table.
Rotor number,N
Revolving speed,ω (rad/s−1)
Equatorial moment of inertia,J×108 (m4)
Unit number,N
Quality,m (kg)
Rotational inertia,Ip(kg⋅m2)
Spring stiffness,k×10−7 (N⋅m−1)
Length,l (m)
1
994.88
3.9760
1
0.325
0
2.5
0.075
2
10.45
0.1125
0.075
3
0.935
0
0.225
4
1.3205
0
0.225
5
0.925
0
0.075
6
8.45
0.09
0.075
7
0.305
0
1.5
0
2
1361.36
5.2002
8
0.1932
0
1.2
0.075
9
7.7796
0.0735
0.15
10
0.3278
0
0.15
11
6.2796
0.0432
0.075
12
0.1732
0
1.0
0In order to establish the transfer matrix of the dual-rotor system, the rotor model in Figure4 is divided into elements as shown in Figure 5, where 10, 11, 17, and 18 are imaginary elements which have no length and mass.Figure 5
Cell partition diagram.The critical speed results of the dual-rotor system calculated by the method in this paper are listed in Table2. The exact solutions are obtained by Lagrange’s equation of motion.Table 2
First six-order critical speed results of various algorithms.
Order
Proposed
Reference
Theoretical
1
400.75
400.97
400.72
2
894.42
894.42
894.41
3
1433.01
1434.26
1432.09
4
1829.76
1829.76
1829.74
5
2119.37
2119.42
2119.30
6
3057.71
3063.46
3056.65
7
3423.41
3429.56
3422.19
8
3723.98
3730.15
3721.94
9
4755.99
4763.05
4752.94
10
12046.95
12077.82
12033.93It can be seen from the results in Table2 that when calculating the first five-order critical speeds, the results of the two methods are close to the exact solutions. However, when exceeding the sixth critical speed, the numerical errors of the algorithm proposed in the literature are obviously increased, while the transfer symplectic matrix method proposed in this paper retains high accuracy.Through the above analysis, it can be found that the difference between the method used in the literature and the method proposed in this paper is very little, and it is difficult to find the advantages and disadvantages of the two methods from the results. Therefore, the relative error analysis is made as follows. The relevant results are shown in Figure6.(a)
From the figure, we can find that the error curves of the two methods show an upward trend overall. The results are consistent with our cognition: when the rotation speed is high, the values of matrix elements may vary greatly, which directly leads to the reduction of calculation accuracy and the increase of error.(b)
When the rotor speed reaches the third-, sixth-, and ninth-order critical speed, the relative errors between the numerical and exact solutions obtained by the two methods are significantly increased. This phenomenon may be caused by the structure of rotor itself.(c)
When the rotor rotates at a low speed (below the fourth critical speed), both the reference method and the method proposed in this paper have satisfactory numerical results in terms of accuracy. However, when the rotor speed exceeds the fifth critical speed, the numerical stability of the method proposed in this paper is obviously better than that adopted in the literature. The relative errors are 0.0347% and 0.2228, respectively, at the sixth critical speed.Figure 6
Comparison diagram of relative error between numerical solution and exact solution.The calculated natural modes are shown in Figure7 (positive vortex).Figure 7
First six formation diagrams. (a)ω1=400.75rad/s. (b) ω2=894.42rad/s. (c) ω3=1433.01rad/s. (d) ω4=1829.76rad/s. (e) ω5=2119.37rad/s. (f) ω6=3057.71rad/s.
(a)
(b)
(c)
(d)
(e)
(f)


## 5. Conclusion

In this paper, the idea of finite element is included into the transfer matrix. Starting from the variational principle, dual variables are introduced to deduce the transfer symplectic matrix of gyro rotor under the Hamilton system.(a)
Compared with the traditional transfer matrix method, a series of complex operations such as solving the aggregation model of beam are avoided, which is easy to understand.(b)
This method is applicable to all kinds of boundary conditions and supporting conditions, so it has universal significance. The state vector reflects the mechanical characteristics of bearing rotor, and the physical meaning is clear and intuitive.(c)
The simulation experimental results show that, in the calculation of high-order critical speed (exceed the fifth critical speed), the numerical accuracy of the transfer symplectic matrix method is obviously better than that of the reference method. Meanwhile, the numerical stability is improved compared with the reference method to some extent, which is reflected by the trend of the curve. The numerical example indicates the feasibility and superiority of the method, which provides the basis for the optimal design of the rotor system.



---

*Source: 9154272-2019-09-09.xml*