Research on modeling and dynamic characteristics of complex coaxial rotor system

To make up deficiency of the finite element method in predicting nonlinear dynamic characteristics of coaxial rotor systems, nonlinear dynamic model of a coaxial rotor system was established with a method combining the finite element method and the fixed interface modal synthesis method. Then an implicit time domain method was presented to solve the nonlinear equations of motion thus dynamic characteristics of the rotor system can be obtained. The computational efficiency of this method largely depends on the number of degrees of freedom with nonlinear forces acting on. With nonlinear forces of squeeze film damper and intermediate bearing considered, nonlinear dynamic response characteristics of the coaxial rotor system under multiple unbalance forces were studied in this work. The results showed that the unbalance excitation frequencies are dominant in the responses of the rotor system. Besides, due to coupling effect of the intermediate bearing some combinations of the unbalance excitation frequencies were also observed in the spectrogram. Stability and periodicity of the rotor system was investigated with bifurcation diagram, Poincare map and phase diagram. It was found that the rotor system executes multiple periods orbital motion under relatively low rotational speeds. With the increasing of rotational speed, the rotor system would execute quasi-periodic motion, chaotic motion and periodic motion again. The quasi-periodic motion and chaotic motion are closely related with the SFD. Finally, under relatively low speed, the nonlinear model was validated by comparing the simulation results with the experimental data. The proposed modeling and solving method is expected to provide theoretical and engineering basis for improving prediction of nonlinear dynamic characteristics of complex rotor systems.


Introduction
Rotor system is the main source of vibration of aeroengines.Stability and reliability of aeroengines mostly depend on dynamic behaviors of the rotor system.In recent years, counter-rotating coaxial system has been applied to many aeroengines, such as RB211, GE120, F119 and so on.The counter-rotating technology is beneficial to improve the fuel consumption rate and thrust-weight ratio as well as reduce the gyroscopic torque of aeroengines [1].Although advantageous for many aspects, it has made dynamic behaviors of the rotor system more complex, especially when coupled with squeeze film dampers(SFD) and intermediate bearing, both of which have been widely applied in modern aeroengines.Surveys of the dynamic issues associated with counter-rotating dual-rotor system are listed in references [2,3].
The most commonly recurring problems in rotor dynamics are excessive steady state synchronous vibration levels and subsynchronous rotor instabilities.The first problem is caused mostly by the synchronous unbalance force and, generally, the second one is caused by nonlinear forces exerted by nonlinear components such as SFD and bearing.The steady state synchronous vibration levels may be reduced by improved balancing, or by moving critical speeds of the system out of the operating range [4], or by introducing external damping to limit peak amplitudes at critical speeds [5].Subsynchronous rotor instabilities can be avoided by rising the natural frequency of the rotor system, or by introducing damping to increase the onset rotor speed of instability [6].However, all these problem solving measures are based on accurate approximation of the actual rotor system.
To better approximate the rotor system, all the forces, linear and nonlinear, should be considered in the dynamic model.However, the problem can be quite complicated when nonlinear forces are considered.The problem of nonlinear response of rotor systems are discussed in the extensive investigations of references [7][8][9][10][11][12][13][14][15].Glasgow D. A. [8] employed the fixed interface modal synthesis method to reduce the size of the rotor system to obtain better computational efficiency.In reference [11,12], nonlinear response of a dual-rotor system with intermediate bearing was studied with transfer matrix method.Wang W. [14] improved the free interface modal synthesis method to analyze dynamic characteristics of a multi-rotor system.In 2004, Hsiao-Wei [15] established the mathematical model of a dual-rotor system with the finite element method to investigate the influences of speed ratio of the rotor system and stiffness of the intermediate bearing.From all the above references it can be seen that the transfer matrix method, finite element method and component modal synthesis method are commonly used in rotor dynamic analysis.In present work, nonlinear model of a counter-rotating coaxial test rig is established with a method combining the finite element method and the component modal synthesis method for accuracy and efficiency.
Notwithstanding the many investigations that there have been of the dynamic behavior of rotor systems with nonlinear forces acting on, the problem of developing efficient method for analyzing nonlinear dynamic behaviors of rotor system still needs to be solved.All these facts necessitate that more detailed models should be established and more efficient and robust numerical methods should be developed to give better approximation of dynamic behaviors of rotor systems.
In this paper, a counter-rotating coaxial system was modeled with finite element method and fixed interface modal synthesis method.Nonlinear forces exerted by the squeeze film damper and the intermediate bearing were considered.With this nonlinear model and the improved Newmark-method, dynamic characteristics of the rotor system was obtained and then further analysis was made and compared with the experimental results.

Finite element formulation
Fig. 1 shows the rotor element considered in this work.Each element consists of two nodes, with four degrees of freedom at each node, namely the lateral displacement and the slope.The lateral displacements in and directions are denoted with and respectively.The slopes are denoted with and for angular displacements in and directions respectively.
The nodal displacement vector can be described as: = . (1) The displacement vector within the element can be interpolated as: The shape functions in Eq. ( 2) are:  In Eq. ( 3), is the element length and is formulated as given below: where the Young's Modulus, the area moment of inertia, the shear modulus, the cross-sectional area. is the shear coefficient, and for a hollow circular shaft section [16,17]: where is the Poisson's ratio, is the ratio of the inner shaft radius to the outer shaft radius.Hence for a solid shaft → 0.
With rotary inertia and shear effects considered, the kinetic energy and strain energy for a shaft element is: where denotes the density of the material and is the transverse shear strain [18].And thus the element mass matrix: The element stiffness matrix: where: The element gyroscopic matrix: where: For disk element, the mass and gyroscopic matrices are: where , and are the mass, diametral moment of inertia and polar moment of inertia of the disk, respectively.

Nonlinear forces of the supports
For squeeze film dampers analyzed in this work, based on the short bearing assumption and the Reynolds boundary conditions, nonlinear forces of the SFD can be expressed as [19,20]: where and are the horizontal and vertical displacements of the journal.( = 1, 2, 3) are Sommerfeld integrals.and are radius and length of the SFD respectively.and denote the dynamic viscosity of the film and the radial clearance of the SFD.
As for the rolling ball bearing, based on pure rolling assumption and the Hertz contact theory, nonlinear force of the rolling ball bearing [21] is: where the superscript and in Eq.(14c) denote inner ring and outer ring of the bearing. in Eq.(14a) is the contact stiffness between rollers and rings.
represent number of the rollers.= 3/2 for the rolling ball bearing used in the intermediate support. is the rotation angle of the th roller at time .and are bearing radial clearance and rotational speed of the bearing retainer respectively.
is the elastic radial deformation of the th roller.and represent radius of the inner and outer ring.and are the rotational speeds of the inner and outer rings.

Equations of motion and numerical algorithm
Equation of motion of a nonlinear rotor system can be written as: where is the unbalance force vector.
is the nonlinear forces exerted by the SFDs and rolling bearing in this work., , and can be obtained quite readily with elements formulated in Section 2.1.Mathematically, the model is a set of nonlinear second order differential equations.The computational efficiency of which depends on the numerical methods applied and the total DOFs of the rotor system.To combine accuracy of the finite element model and computational efficiency, the fixed interface modal synthesis is applied to reduce dimension of the mathematical model and thus the computational effort.
Eq. ( 15) can be rewritten as: where can be interpreted as linear DOFs, with only unbalance forces acting on, while represents the interface DOFs or nonlinear DOFs in this work, including all the nodes with nonlinear forces acting on.
According to fixed interface modal synthesis method, transformation between physical and modal coordinates is: where -the mass normalized normal mode matrix, -the mass normalized constrained mode, -unity matrix, -the normal mode coordinates, -the transformation matrix. is given as below: In Eq. ( 16), let = 0 and neglecting the gyroscopic effects gives: Solving the eigenvalue problem corresponding to Eq. (17c) gives the mass normalized normal mode matrix and the modal frequency Ω .Substitute Eq. (17a) into Eq.( 16) and premultiply with to obtain the reduced equations of motion: where: where Ω is the th modal frequency obtained from Eq. (17c). is the unbalance force acting on linear DOFs. is the nonlinear force acting on nonlinear DOFs.And: where = ⋯ is the nonlinear force vector exerted by SFDs which can be calculated with Eq. ( 13).= − − is the nonlinear force vector exerted by rolling bearing which can be calculated with Eq. (14).and are given below: where and are the stiffness and damping coefficients of the elastic support, respectively.is the number of elastic supports.
Rearranging Eq. ( 18) gives: In Eq. ( 23), the vector and can be interpreted as linear and nonlinear DOFs of the reduced system respectively.Obviously, there is no nonlinear force acting on the linear DOFs.Thus, the explicit Newmark-method applies to Eq.(23a) while implicit Newmark-method applies to Eq. (23b).
According to assumptions of Newmark-method, in time interval : where = + Δ and Δ is the time increment.The superscript and + 1 denote and .Following equations can be obtained from Eq. (24): With: Substituting Eq. ( 24)-Eq.( 27) into Eq.( 23) yields: With: Using Eq. ( 26) to substitute for in Eq. (29g) and substitute the new Eq.(29g) into Eq.(28b) to yield nonlinear equations of .After solving Eq. (28b) for with numerical algorithms, can be obtained with substituted into Eq.(28a).Thus , , , and can be solved from Eq. (25) and Eq. ( 26).The process continues to repeat and move forward to find , and so on.When the iteration is over, the nonlinear response in physical coordinate system can be obtained from Eq. (17a).
Computational efficiency of solving Eq. ( 23) depends on solving Eq. (23b) while dimensions of Eq. (23b) depend on DOFs with nonlinear forces acting on.Thus, the computational efficiency can be greatly improved with the modeling technique and solving method described in this work.
To summarize, nonlinear model of the rotor system is established with finite element method and fixed interface modal synthesis method; subsequently an implicit time-domain method based on Newmark-method is applied to solve the equations of motion of the reduced system thus dynamic characteristics can be obtained.The process of modeling and solving is described as follows: First, only the rotating components of the rotor assembly are modeled with the element described in Section 2.1 thus , and in Eq. ( 15) can be obtained and partitioned as Eq. ( 16) for calculating and later use.Then, in accordance with Eq. (17c), modal analysis is conducted on the model with all supports pinned to obtain and Ω .The left hand side of Eq. (18a) now is obtained with , and the partitioned , and while the right hand side is formulated with Eq. ( 13), Eq. ( 14) and Eq. ( 19)-Eq.( 22).Finally, Eq. (18a) is solved with the numerical method which is described by Eq. (23)-Eq.(29).Flowchart of the modeling and solving method is shown in Fig. 2.

Test rig description
An overall view of the dual-rotor test rig is given in Fig. 3.The test rig was designed to study the dynamic characteristics of co-and counter-rotating dual-rotor systems with 4 or 5 supports.In this work, the 4 supports counter-rotating dual-rotor system is studied.So, the additional support in Fig. 3 was removed.Thus, the studied rotor system is designed with 4 supports and 4 disks, two for each rotor.The test rig consists of two shafts disposed along the same axis, connected by an intermediate bearing.Each shaft is driven by its own motor thus their spin speeds can be different.Speed ratio of outer and inner rotor is -1.65.Rear support of the outer rotor is intermediate support in which a rolling ball bearing is adopted.Support I and II are front and rear supports of inner rotor respectively while Support III and IV are front and rear support of outer rotor.The squirrel cage + rolling bearing + SFD supporting scheme is adopted for support I, II and III.While support IV, or the intermediate bearing, is mounted without SFD.Four eddy current probes are placed to measure the displacement of disk 2 and 4, two for each disk.It is to be noted that in Fig. 3 the eddy current probes have not been placed in the right positions yet.
The test profiles include slow run-up and run-down through the speed range of the test rig.The angular acceleration for inner rotor during run-up and rundown is approximately 10 rpm/s.Besides, the steady state response are measured when inner rotor spins at 128 rad/s, 160 rad/s and 210 rad/s.The operational ranges for the inner and outer rotor are 0-232 rad/s and 0-382 rad/s respectively.The excitation is only due to the residual unbalance of the rotors.

Model introduction
The counter-rotating coaxial system studied in this paper is presented in Fig. 4. The stationary coordinate system in Fig. 4 consists of three mutually perpendicular axes, , , , intersecting at the point and axis of the rotor coincides with axis .The axes , are horizontal, is vertical.Model parameters of the rotor system studied in this work are listed in Table 1-5.The Young's modulus of the shaft is 196 GPa, mass density 7810 kg/m 3 , and shear modulus 75.5 GPa.To apply modal synthesis method, 40 modes are retained in the normal mode matrix .For the Newmark-method, = 0.25, = 0.5.
Geometric dimensions and information of each element are listed in Table 1.Stiffness coefficients of elastic supports are listed in Table 2. Parameters of the intermediate bearing and SFDs are listed in Table 3 and 4

Unbalance response analysis
With the model established by method described in Section 2, steady-state unbalance response of the rotor system is obtained for rotational speed of the inner rotor varying from 4 rad/s to 1600 rad/s with a step length of 2 rad/s.Zero initial condition is adopted for the first step.For the rest steps, result of previous step is adopted as initial condition.Without loss of generality, unbalance response of disk 2 and disk 4 are analyzed to evaluate the coupling effect of the intermediate bearing.Spectrum cascades of the horizontal response of disk 2 and disk 4 with inner rotor operating in 4-400 rad/s are shown in Fig. 5 and Fig. 6.The simulation and experimental results are shown in Fig. 5 and Fig. 6 respectively.The numerically simulated spectrum cascades of the horizontal response of disk 2 and disk 4 with inner rotor operating in 400-1600 rad/s are shown in Figs.7-9.Due to nonlinearities of SFD and the intermediate bearing as well as the multifrequency unbalance excitation of the inner and outer rotor, rich frequency components are observed in the response of disk 2 and disk 4. Following conclusions can be reached from Fig. 5-9: (1) Responses of the inner and outer rotor are coupled because of the intermediate bearing.Frequency components, and , corresponding to the unbalance excitation frequencies of the inner and outer rotor are dominant in the response of disk 2 and disk 4. The coupling effect also makes the spectrum cascades corresponding to disk 2 and disk 4 basically follow the same pattern with only slight differences in amplitude.
(3) By comparison between Fig. 5 and Fig. 7 it can be seen that the new combination frequency components, 2 − , 3 − 2 emerge when inner rotor operating in 400-800 rad/s.However, 2 − and 3 − 2 are not observed for both disk 2 and 4 in Fig. 8 and Fig. 9 when the inner rotor operating in 800-1200 rad/s and 1200-1600 rad/s.Besides, 2 + is persistent in 0-1600 rad/s for disk 2 and 4, which can be seen from Fig. 5 and Figs.7-9.
(4) Two critical speeds of the rotor system, 168 rad/s and 186 rad/s, are observed in rotational speed range 0-200 rad/s in Fig. 5.The first one, 168 rad/s, is excited by the outer rotor and the second one, 186 rad/s, is excited by the inner rotor.It can be obtained from the experimental data shown in Fig. 6 that the corresponding critical speeds of the rotor system are 173 rad/s and 189 rad/s.Discrepancies of critical speeds between numerical and experimental results are less than 3 %, which demonstrate the validity of the model established in this work.Besides the pattern of transformation between different states of motion, it is found that the chaotic motion is closely related with the SFD.Fig. 11(a), Fig. 15(a) and Fig. 18(a) show the spectrogram of chaotic motion at different rotational speeds, the same low frequency component, approximately 30 Hz, can be seen quite obviously.This frequency component, 30 Hz or 188 rad/s, is almost the value of the first critical speed excited by the inner rotor.This phenomenon suggests that the SFD has an important influence on the state of motion of the rotor system.

Experimental validation
The comparison between numerical and experimental results has already been made in previous section, Fig. 5 and 6.Orbits of disk 2 and 4 under three different rotational speeds are compared in this section for further validation and analysis, see Fig. 19  Following conclusions can be reached by analysis of Fig. 19-Fig.21: (1) With dual unbalance excitation, disk 2 and 4 execute petal-shaped orbits.Especially when approaching critical speeds of the rotor system, Fig. 20.
(2) The rotor system is slightly anisotropic, which might be caused by the squirrel cage, see Fig. (3) Under relatively low speed, the numerical results show great agreement with the experimental results.

Conclusions
A modeling method combining the finite element method and the fixed interface modal synthesis method has been developed in the work to establish the nonlinear model of a counter-rotating dual-rotor test rig for steady state nonlinear response analysis.First, a finite element model of the rotor system is established without considering supports, from which the mass, stiffness and gyroscopic matrices are obtained.Together with these matrices, the fixed interface modal synthesis method is applied to establish the nonlinear model of the rotor system in which the nonlinearities of SFD and intermediate bearing are considered.Subsequently, the Newmark-method is improved to solve the nonlinear equations of motion.Then dynamic characteristics of the rotor system are investigated.
Conclusions are listed as following: 1) The modeling method developed in this work is fast and efficient in establishing nonlinear model of complex rotor systems.Combing with the improved Newmark-method, nonlinear dynamic characteristics can be obtained accurately and efficiently.
2) Due to coupling effect of the intermediate bearing, dual unbalance excitation frequencies are dominant in the responses of inner and outer rotor.Besides, frequency components other than dual unbalance excitation frequencies in the responses of inner and outer rotor are mainly caused by coupling effect of the intermediate bearing.
3) Rotors execute stable periodic orbital motion in 4-390 rad/s.Transformations between periodic, quasi periodic and chaotic motion are reached with increasing rotational speed.A low frequency component which is approximately the first critical speed of the outer rotor is observed, which means the SFD has a great influence on the state of motion of the rotor system.4) Orbit of the rotor system under dual unbalance excitation is petal shaped, especially when approaching critical speeds of the rotor system.Under relatively low speed, the validities of the modeling method and the solving method were demonstrated with the experimental results.

Fig. 4 .
Fig. 4. Structural diagram of the coaxial rotor system

a) Disk 2 b) Disk 4 Fig. 6 .
Spectrum cascade of the horizontal response of disk 2 and 4-experimental results

a) Disk 2 b) Disk 4 Fig. 7 .b) Disk 4 Fig. 8 .b) Disk 4 Fig. 9 .
Spectrum cascade of the horizontal response of disk 2 and 4-numerical simulation a) Disk 2 Spectrum cascade of the horizontal response of disk 2 and 4-numerical simulation a) Disk 2 Spectrum cascade of the horizontal response of disk 2 and 4-numerical simulation

Fig. 10 4 Fig. 10 . 11 . 12 . 13 . 15 . 16 . 17 .
Fig. 10 shows bifurcation diagrams for horizontal displacement of disk 2 and disk 4 with rotor speed.Sampling period for the bifurcation diagram is 5×2 / because the rotational speed ratio of outer and inner rotor is -1.65 and 5 ≈ 1.65×3.The horizontal axes of Fig. 10(a) and (b) are rotational speeds of inner and outer rotor respectively.And the vertical axes are horizontal displacement of disk 2 and disk 4.Basically, disk 2 and 4 execute multiple period orbital motion with inner rotor operating in 4-390 rad/s, see Fig.10(a), Fig.10(b).However, with the increase of rotational speed, a transformation between periodic and chaotic motion can be observed.From 400 rad/s to 600 rad/s, .
a) Disk 2: numerical results b) Disk 2: experimental results c) Disk 4: numerical results d) Disk 4: experimental results

Table 1 .
. The unbalance configuration and inertia properties of each disk are listed in Table5.Dimension and elements information of the rotor system Node no.Axial location (m) Bearing/disk Element no.

Table 3 .
Parameters of the intermediate bearing

Table 5 .
Unbalance configuration and inertia properties of disks