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

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

Figure 1 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 plane of the discrete rotor model, a typical unit consists of a massless shaft segment and a lumped mass disc [19].

**(a)**

**(b)**

Considering the element, defining as the stator vector containing deflection , with as the slope, as the bending moment, and as the shear force, we obtain

For any axial segment , there is a certain relationship between the state vectors , that iswhere 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 node represent the bending moment and shear force to the left side of axis segment , respectively, while the bending moment and shear force to the left of node represent the bending moment and shear force to the right side of axis segment , respectively.

According to the bending deformation formula and displacement analysis of beam, the following equation is obtained [20–22]:where the field transfer matrix isin which represents the shear effect coefficient, is called the cross-sectional shape coefficient, and 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 . Suppose that the node has a spring bearing and the spring stiffness is . 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]where the transfer matrix isin which is the moment of inertia of diameter, 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 from to can be written as

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):where the assembled transferring matrix is

It is found that the value of the element containing in the total transfer matrix 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 vector is converted intowhere and represent the force part of the state vector and the displacement part of the state vector, respectively. A generalized Riccati transformation at the section is given as follows:where is the so-called Riccati transfer matrix. And through the derivation, we get the recursion formula

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

Consider a dual-rotor system with two intermediate support, as shown in Figure 3. Defining the state vectors of the element of the two rotors as , respectively:

As for the noncoupling element, according to the transfer matrix method of single rotor system, their respective transfer relations can be obtained:denoted byin which are the transfer matrix of the two rotors. Their form can be derived from the Prohl, Riccati, or other transfer matrix methods. , , and , where are defined as 0 when there is no spring bearing in the element.

As for the coupling element, assume that the coupling point is in the unit. According to the balance of forces, the coordination equations are obtained:so that we can get state vectors on both sides of the coupling element:where represents the transfer matrix for the coupling element, denoted aswhere stands for support stiffness of the spring and is the torque stiffness.

##### 2.3. Critical Speed and Mode of Vibration

Supposing that the assembled transferring matrix is denoted by , according to the boundary conditions of the model

Combine equations (17) and (19), there is

Equation (20) is a system of homogeneous algebraic equations, where each element is a function of vortex frequency, and its coefficient determinant called surplus is

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:

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 the *z*-coordinate be along the axis of the beam before deformation. Let be the cross-section shear coefficient, be the shear modulus, be the elastic modulus, be the moment of inertia, be the cross-section area of the beam, and and be the linear and angular displacements of the cross section, respectively. The deformation energy of this beam can be expressed as

Take the perpendicular plane, for example, the internal force of the beam can be expressed as

The dynamic equations of the beam can be deduced to

For the sake of analysis, the time domain is converted to the frequency domain and the variable is still represented by the previous notationwhere is the rotation angular velocity of the rotor, and then, the dynamic equation (25) can be converted into

By the same logic, we can get that of the plane. Thus, we can obtain the kinetic energy of this beam:

Therefore, the Lagrangian density function is

For the convenience of analysis and calculation, the following notation is introduced: . Take the matrix form

Then, the potential energy variation principle can be given for the above equation as follows:

The corresponding variational principle is and of which

Then, the dynamic equation can be re-expressed as

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, let and then it can be re-expressed as

Substitute equation (34) into the dynamic equation, we can getwhere .

Then, the Hamiltonian is introduced as

The corresponding variational principle is , and the dual vectors are unified to form the full state vector .

The equation above is for the general -dimensional displacement vector; in this paper, there is , then . Moreover, at the same time, we can get the Hamiltonian equation of state asin which

It is observed that the mechanical meaning of dual vector is composed of generalized internal force, namely, shear force and bending moment , which shows the rationality and validity of the introduced dual variables.

When solving the ordinary differential equation (37), we let and

Along the direction of the -axis, divide the -axis into arbitrarily small intervals. The transformation matrix when transferring from section to section is

Now that the is small enough, we can expand by Taylor series

Equation (41) indicates that the series converges so that the method put forward is available. Then, , in which is the transfer matrix 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 Figure 4. Some material and geometrical properties are given as *E* = 206 GPa, *μ* = 0.3, and *ρ* = 7850 kg/m^{3}. The detailed original data of the calculation example are shown in Table 1.

In order to establish the transfer matrix of the dual-rotor system, the rotor model in Figure 4 is divided into elements as shown in Figure 5, where 10, 11, 17, and 18 are imaginary elements which have no length and mass.

The critical speed results of the dual-rotor system calculated by the method in this paper are listed in Table 2. The exact solutions are obtained by Lagrange’s equation of motion.

It can be seen from the results in Table 2 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 Figure 6.(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.

The calculated natural modes are shown in Figure 7 (positive vortex).

**(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.

#### Data Availability

The data used to support the findings of this study have not been made available because the data that have been used are confidential and the authors do not have permission to share data.

#### Conflicts of Interest

The authors declare no potential conflicts of interest with respect to the research, authorship, and/or publication of this article.

#### Acknowledgments

This work was supported by the National Natural Science Foundation of China (Grant nos. 61573012 and 61702388) and Fundamental Research Funds for Central University (WUT2016A004).