Abstract
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. With nonlinear forces of squeeze film damper and intermediate bearing considered, nonlinear dynamic response characteristics of the co- and counterrotating coaxial rotor system under multiple unbalance forces were studied and compared in this work. It was found that the critical speeds of the corotating system were equal to or slightly higher than those of the counterrotating case due to the gyroscopic moments. 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 but the combinations were different for co- and counterrotating cases. Stability and periodicity of the rotor system were investigated with bifurcation diagram, Poincare map, and phase diagram. It was found that the rotor system executes four-period quasi-periodic motion around critical speeds.
1. Introduction
In recent years, coaxial dual-rotor system has been applied to many aeroengines with two corotating or counterrotating rotors, such as RB211, GE120, and F119. The counterrotating technology is beneficial to improve the fuel consumption rate and thrust-weight ratio as well as reduce the gyroscopic torque of aeroengines [1]. Another technology that has been applied is the intermediate bearing which aims to further reduce the mass of the aeroengines. For most modern aeroengines, the squeeze film dampers (SFD) are commonly adopted to provide structural isolation and reduce the amplitude of rotor response. All these facts necessitate that more studies should be conducted to investigate the dynamic characteristics of dual-rotor systems with SFD and the intermediate bearing.
In 1975, Vance and Royal [2] have published an extensive discussion of the design and operational technological issues related to the intermediate bearings. Hibner [3] has applied the transfer matrix method to the multiple-shaft machines in order to compute the critical speeds and nonlinearly damped response. Gupta et al. [4] have presented a study on a counterrotating dual-rotor test rig with an intermediate bearing. The rotor system was also modeled with the transfer matrix method. Cross-excitation phenomena have been encountered. In 1996, Ferraris et al. [5] have analyzed the rotordynamics of a dual-shaft prop-fan aircraft engine with two rotors spin at equal speeds in opposite directions. A study of a twin-spool aircraft engine is also presented in the book of Lalanne and Ferraris [6]. In [7–9], Hu et al. have studied the dynamic characteristics of a counterrotating dual-rotor system with the transfer matrix method. The nonlinearity of the intermediate bearing was included. The finite element method has been applied in [10, 11] to obtain the dynamic characteristics of dual-rotor system with intermediate bearing. Fatigue life of two different configurations of the intermediate bearing were studied and compared by Hu et al. [12] in 2006. However, the intermediate bearing was not coupled with the rotor system. From 2008 to 2011, Hai and Bonello [13–16] have conducted extensive investigations on computational method to obtain the unbalance response of the multispool whole aeroengine with squeeze film dampers.
The current study is on a co- and counterrotating dual-rotor test rig that has been developed to study the dynamics of dual-rotor machines with intermediate bearings and squeeze film dampers. The purpose of this paper is to study the modeling method and nonlinear dynamic characteristics of the complex dual-rotor systems. In the second part of the paper, the modeling method and the numerical algorithm are presented. Then a description of the experimental apparatus is given in the next part. Finally, nonlinear dynamic characteristics of the co- and counterrotating test rig are studied and the numerical results are compared with the experimental results to validate the model established in this paper.
2. Modeling and Numerical Algorithm
2.1. Finite Element Formulation
Figure 1 shows the rotor element considered in this work. Each element consists of two nodes, with four degrees of freedom at each node. The nodal displacement vector can be described asThe displacement vector within the element can be interpolated asThe shape functions in (2) areIn (3a) and (3b), is the element length and is formulated as given below:where is Young’s Modulus, is the area moment of inertia, is the shear modulus, and is the cross-sectional area. denotes the shear coefficient, and for a hollow circular shaft section [16, 17]where is Poisson’s ratio and is the ratio of the inner shaft radius to the outer shaft radius. Hence for a solid shaft .
With rotary inertia and shear effects considered, the kinetic energy and strain energy for a shaft element areAnd thus the element mass matrix isThe element stiffness matrix iswhereThe element gyroscopic matrix iswhereFor disk element, the mass and gyroscopic matrices arewhere , , and are the mass, diametral moment of inertia, and polar moment of inertia of the disk, respectively.
2.2. 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 [17] and are the horizontal and vertical displacements of the journal. 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 [18] iswhere the superscripts ir and or in (14c) denote inner ring and outer ring of the bearing. in (14a) is the contact stiffness between rollers and rings. represent number of the rollers. for the rolling ball bearing used in the intermediate support. is the rotation angle of theth 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.
2.3. Equations of Motion and Numerical Algorithm
Equation of motion of a nonlinear rotor system can be written as 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.
Equation (15) can be rewritten aswhere can be interpreted as linear DOFs, with only unbalance forces considered, while represents the interface DOFs or nonlinear DOFs in this work, including all the nodes with nonlinear forces considered.
According to fixed interface modal synthesis method, transformation between physical and modal coordinates iswhere is the mass normalized normal mode matrix, is the mass normalized constrained mode, is unity matrix, is the normal mode coordinates, and is the transformation matrix. is given as below:In (16), let and neglecting the gyroscopic effects givesSolving the eigenvalue problem corresponding to (17c) gives the mass normalized normal mode matrix and the modal frequency .
Substitute (17a) into (16) and premultiply with to obtain the reduced equations of motion:where is the th modal frequency obtained from (17c). is the unbalance force acting on linear DOFs. is the nonlinear force acting on nonlinear DOFs.
Andwhere is the nonlinear force vector exerted by SFDs which can be calculated with (13a), (13b), (13c), (13d), and (13e). is the nonlinear force vector exerted by rolling bearing which can be calculated with (14a), (14b), (14c), (14d), and (14e). and are given below: and are the stiffness and damping coefficients of the elastic support, respectively. is the number of elastic supports.
Rearranging (18a) and (18b) givesIn (21a) and (21b), the vectors 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-beta method applies to (21a) while implicit Newmark-beta method applies to (21b).
According to assumptions of Newmark-beta method, in time interval iswhere and is the time increment. The superscripts and denote and .
The following equations can be obtained from (22a) and (22b):withSubstituting (22a) and (22b)–(25) into (21a) and (21b) yieldswithUse (24) to substitute for in (27g) and substitute the new (27g) into (26b) to yield nonlinear equations of . After solving (26b) for with numerical algorithms, can be obtained with substituted into (26a). Thus , , , and can be solved from (23) and (24). 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 (17a).
Computational efficiency of solving (21a) and (21b) depends on solving (21b) while dimensions of (21b) depend on DOFs with nonlinear forces considered. 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-beta method is applied to solve the equations of motion of the reduced system; thus dynamic characteristics can be obtained. Flowchart of the modeling and solving method is shown in Figure 2.
3. Test Rig Description
The coaxial test rig studied in this paper is presented in Figure 3.
The stationary coordinate system in Figure 3 consists of three mutually perpendicular axes, , , and , intersecting at the point and axis of the rotor coincides with axis . The axes and are horizontal, and is vertical.
The rotor system consists of two shafts disposed along the same axis , connected by an intermediate bearing. The test rig is designed with 4 supports and 4 disks, two for each rotor. The rotors are supported by supports I, II, III, and IV. The squirrel cage + rolling bearing + SFD supporting scheme is adopted for supports I, II, and III. The intermediate bearing has one ring mounted on the shaft of inner rotor and the other ring on outer rotor as shown in Figure 3. The rotor comprises shaft and disks. Each shaft is driven by its own motor through flexible couplings and therefore their rotational speeds can be different.
Model parameters of the rotor system studied in this work are listed in Tables 1–5. Young’s modulus of the shaft is 196 GPa, mass density is 7810 kg/m3, and shear modulus is 75.5 GPa. To apply modal synthesis method, 40 modes are retained in the normal mode matrix . For the Newmark-beta method, and .
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 Tables 3 and 4. The unbalance configuration and inertia properties of each disk are listed in Table 5.
4. Numerical and Experimental Results
Dynamic characteristics of both the co- and counterrotating dual-rotor system were studied in this section. The rotational speed ratio is defined as , where is the rotational speed of the outer rotor and is the rotational speed of the inner rotor. (corotating) and (counterrotating) cases are studied in the current work.
4.1. Overview of the Results
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 800 rad/s with a step length of 2 rad/s. Zero initial condition is adopted for the first step. For the rest of the steps, result of previous step is adopted as the initial condition. The critical speeds of the rotor system can be identified by plotting and analyzing the spectrum cascade of the numerical and experimental results. Without loss of generality, only unbalance response of disk 2 and disk 4 are analyzed.
Numerical results of the spectrum cascades for the horizontal response of disk 2 and disk 4 with inner rotor operating in 4–400 rad/s are shown in Figures 4 and 5. And the corresponding experimental results are shown in Figures 6 and 7. The following conclusions can be reached from Figures 4–7:(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 which makes the spectrum cascades corresponding to disk 2 and disk 4 basically follows the same pattern with only slight differences in amplitude.(2)Frequency components besides and in the responses of disks 2 and 4, such as , , , , and , are mainly caused by the nonlinear forces of the SFD and the coupling effect of the intermediate bearing. However, the combination frequency components are different for co- and counterrotating systems. For example, Figures 4 and 5 show the spectrum cascade of disks 2 and 4 under counter- and corotating conditions. In Figure 4(a), the combination frequency components are , , , and , while in Figure 5(a) only and are observed. In Figure 4(b), the combination frequency components are , , , , and while in Figure 5(b) those combinations are ,, and.(3)In Figure 4, two critical speeds of the counterrotating rotor system, 168 rad/s and 186 rad/s, are observed. For the corotating system, Figure 5, the corresponding critical speeds are 168 rad/s and 190 rad/s. From the experimental results presented in Figures 6 and 7, the lowest two critical speeds for counter- and corotating systems are 170 rad/s and 189 rad/s and 173 rad/s and 195 rad/s. The differences of critical speeds for co- and counterrotating systems are mainly caused by the gyroscopic effect. For the corotating rotor system, the gyroscopic moments of both rotors act in the same direction to strengthen the shaft thus raising the critical speeds. While, for the counterrotating system, the gyroscopic moments act in the opposite directions; thus the strengthening effect is weakened and then lower critical speeds are achieved.
(a) Disk 2
(b) Disk 4
(a) Disk 2
(b) Disk 4
(a) Disk 2
(b) Disk 4
(a) Disk 2
(b) Disk 4
4.2. Unbalance Response Analysis
Respectively, for counter- and corotating cases, Figures 8 and 9 show bifurcation diagrams for horizontal displacement of disk 2 and disk 4 with rotor speed. Sampling period for the bifurcation diagram is because the rotational speed ratio of outer and inner rotor is and 5 ≈ 1.65 × 3. The horizontal axes are rotational speeds of inner and outer rotor. And the vertical axes are horizontal displacement of disk 2 and disk 4.
(a) Disk 2
(b) Disk 4
(c) Disk 2
(d) Disk 4
(a) Disk 2
(b) Disk 4
(c) Disk 2
(d) Disk 4
Basically, disks 2 and 4 execute multiple period orbital motion around critical speeds, which can be seen from the bifurcation diagrams shown in Figures 8 and 9. For both the counter- and corotating cases, when the rotor systems operate near the critical speeds, a four-period orbital motion is observed; see Figures 10, 11, 16, and 17. Although only four scattered and isolated points are observed in the Poincare maps shown in Figures 10(d), 11(d), 16(d), and 17(d), the enlarged graphic circulated by the red circle shows that the isolated “points” are actually closed loops which indicate quasi-periodic motion. Besides, by comparison between Figures 10 and 11 it can be seen that both the counter- and corotating systems execute four-period orbital motion although their orbits show a little difference—the amplitude for the orbit of disk 2 of the counterrotating system is slightly larger than that of the corotating system. Also, only and are observed in the spectrogram shown in Figures 10(a) and 11(a). And the amplitude of is much larger than that of , which means the outer rotor is the main source of vibration.
(a) Spectrogram
(b) Orbit
(c) Phase trajectory
(d) Poincare map
(a) Spectrogram
(b) Orbit
(c) Phase trajectory
(d) Poincare map
As the rotational speed rises to 114 rad/s, notable differences for the spectrogram, orbit, phase trajectory, and Poincare map between the counter- and corotating conditions are observed, Figures 12 and 13. The Poincare map shown in Figure 12(d) suggests that disk 2 executes four-period quasi-periodic motion under counterrotating condition while no obvious periodicity can be observed from the random-like and irregularly distributed points shown in Figure 13(d), the corotating condition. In contrast, this phenomenon is reversed when the rotational speed is further increased to 140 rad/s; see Figures 14 and 15. Under 140 rad/s, the corotating case, Figure 15(d), indicates a four-period quasi-periodic motion while the counterrotating case, Figure 14(d), shows randomly points cloud. Besides, from Figures 13(a) and 15(a) it can be seen that frequency components other than and can be found for the conditions when no obvious periodicity can be observed.
(a) Spectrogram
(b) Orbit
(c) Phase trajectory
(d) Poincare map
(a) Spectrogram
(b) Orbit
(c) Phase trajectory
(d) Poincare map
(a) Spectrogram
(b) Orbit
(c) Phase trajectory
(d) Poincare map
(a) Spectrogram
(b) Orbit
(c) Phase trajectory
(d) Poincare map
(a) Spectrogram
(b) Orbit
(c) Phase trajectory
(d) Poincare map
(a) Spectrogram
(b) Orbit
(c) Phase trajectory
(d) Poincare map
At 170 rad/s, approximately the next critical speed, four-period quasi-periodic motion emerges again; see Figures 16 and 17. Also, it is noted that the closed loop shown in Figure 16(d), the counterrotating case, is twisted while that of the corotating case, Figure 17(d), is an ellipse. Figures 16 and 17 also reveal that the amplitude of disk 2 under counterrotating condition is approximately twice that for the corotating condition.
To conclude, the rotors execute four-period quasi-periodic motion around critical speeds, Figures 10, 11, 16, and 17. Between two critical speeds, the motion state of the rotor system can be quite different for the counter- and corotating cases, Figures 12–15. Around critical speeds, the amplitude of disk 2 is much larger for the counterrotating case comparing with that of the corotating case, Figures 10(b), 11(b), 16(b), and 17(b).
4.3. Experimental Validation
The comparison between numerical and experimental results has already been made in the previous section, Figures 4–7. Orbits of disk 2 and 4 under three different rotational speeds are compared in this section for further validation; see Figures 18–22. The numerical results show great agreement with the experimental results, which validate the model and method adopted in the current work.
(a) Disk 2—numerical results
(b) Disk 2—experimental results
(c) Disk 4—numerical results
(d) Disk 4—experimental results
(a) Disk 2—numerical results
(b) Disk 2—experimental results
(c) Disk 4—numerical results
(d) Disk 4—experimental results
(a) Disk 2—numerical results
(b) Disk 2—experimental results
(c) Disk 4—numerical results
(d) Disk 4—experimental results
(a) Disk 2—numerical results
(b) Disk 2—experimental results
(c) Disk 4—numerical results
(d) Disk 4—experimental results
(a) Disk 2—numerical results
(b) Disk 2—experimental results
(c) Disk 4—numerical results
(d) Disk 4—experimental results
5. 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 counterrotating 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-beta method is improved to solve the nonlinear equations of motion. Then dynamic characteristics of the rotor system are investigated.
Conclusions are listed as follows:(1)The modeling method developed in this work is fast and efficient in establishing nonlinear model of complex rotor systems. Combining with the improved Newmark-beta method, nonlinear dynamic characteristics can be obtained accurately and efficiently.(2)Due to the gyroscopic effect, critical speeds of the corotating rotor system are slightly higher than those of the counterrotating system.(3)Due to coupling effect of the intermediate bearing, the dual unbalance excitation frequencies are dominant in the responses of inner and outer rotor for both the counter- and corotating systems. Besides, combination frequency components in the responses of inner and outer rotor are mainly caused by coupling effect of the intermediate bearing.(4)Generally the rotors execute four period quasi-periodic motion around critical speeds. Between two critical speeds, the motion state of the rotor system can be quite different for the counter- and corotating cases. Around critical speeds, the amplitude of disk 2 is much larger for the counterrotating case comparing with that of the corotating case.(5)The comparison between numerical and experimental results has shown great agreement, which proves validities of the modeling method and the numerical results.
Conflicts of Interest
The authors declare that there are no conflicts of interest regarding the publication of this paper.