# A Comparison Study on Co- and Counterrotating Dual-Rotor System with Squeeze Film Dampers and Intermediate Bearing

Authors: Fei Wang; Gui-Huo Luo; Shu Yan; Hai-Tao Cui
Journal: Shock and Vibration
 (2017)

---

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


---

## Body

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

Figure1 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 as(1)qe=ue1ve1θe1φe1ue2ve2θe2φe2T.The displacement vector within the element can be interpolated as(2)ueξξ,tveξξ,tθeξξ,tφeξξ,t=N100N2N300N40N1-N200N3-N400-Nˇ1Nˇ200-Nˇ3Nˇ40Nˇ100Nˇ2Nˇ300Nˇ4qe.The shape functions in (2) are(3a)N1=1-1le2+12g12gleξ+3ξ2-2leξ3,N2=1le2+12gle2+6gξ-2le2+6gleξ2+ξ3,N3=1le2+12g12gleξ+3ξ2-2leξ3,N4=1le2+12g-6gξ-le2-6gleξ2+ξ3,(3b)Nˇ1=1le2+12g6leξ2-6ξ,Nˇ2=1le2+12gle2+12g-4le2+12gleξ+3ξ2,Nˇ3=1le2+12g6ξ-6leξ2,Nˇ4=1le2+12g12g-2le2leξ+3ξ2.In (3a) and (3b), le is the element length and g is formulated as given below:(4)g=EeIeκGeAe,where Ee is Young’s Modulus, Ie is the area moment of inertia, Ge is the shear modulus, and Ae is the cross-sectional area. κ denotes the shear coefficient, and for a hollow circular shaft section [16, 17](5)κ=61+μ1+λ227+12μ+4μ21+λ22+45+6μ+2μ2λ2,where μ is Poisson’s ratio and λ is the ratio of the inner shaft radius to the outer shaft radius. Hence for a solid shaft λ→0.Figure 1
Diagram of beam element.With rotary inertia and shear effects considered, the kinetic energy and strain energy for a shaft element are(6a)Te=12∫0leρeAeu˙eξ2ξ,t+ρeIeφ˙eξ2ξ,tdξ,(6b)Ue=12∫0le∂φeξξ,t∂ξTEeIe∂φeξξ,t∂ξdξ+12κ∫0leγTGeAeγdξ.And thus the element mass matrix is(7)me=ρeAe∫0leNTNdξ+ρeIe∫0legN′′′+N′TgN′′′+N′dξ.The element stiffness matrix is(8)ke=EeIe∫0leN′′TN′′dξ+EeIeg∫0leN′′′TN′′′dξ,where(9)N′=dN1N2N3N4dξ,N′′=d2N1N2N3N4dξ2,N′′′=d3N1N2N3N4dξ3.The element gyroscopic matrix is(10)ge=-2ρeIe∫0leB2TB1-B1TB2dξ,where(11)B1B2=0-Nˇ1Nˇ200-Nˇ3Nˇ40Nˇ100Nˇ2Nˇ300Nˇ4.For disk element, the mass and gyroscopic matrices are(12a)med=mdi0000mdi0000Idi0000Idi,(12b)ged=00000000000Ipi00-Ipi0,where mdi, Idi, and Ipi 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](13a)fxsfd=-μsRL3c2x2+y21/2xε˙I2+εψ˙I1-yε˙I1+εψ˙I3,fysfd=-μsRL3c2x2+y21/2yε˙I2+εψ˙I1+xε˙I1+εψ˙I3,(13b)ε=x2+y2c,(13c)ε˙=xx˙+yy˙cx2+y2,(13d)ψ˙=yx˙-xy˙x2+y2,(13e)tan⁡ψ=xy.x and y are the horizontal and vertical displacements of the journal. Ijj=1,2,3 are Sommerfeld integrals. R and L are radius and length of the SFD, respectively. μs and c 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] is(14a)fxb=kn∑j=1NbuθjξHuθjsin⁡θj,fyb=kn∑j=1NbuθjξHuθjcos⁡θj,(14b)θj=2πj-1Nb+ωct,(14c)uθj=xir-xorcos⁡θj+yir-yorsin⁡θj-γ2,(14d)Huθj=0,uθj≤0,uθj,uθj>0,(14e)ωc=ωinr+ωoutRR+r,where the superscripts ir and or in (14c) denote inner ring and outer ring of the bearing. kn in (14a) is the contact stiffness between rollers and rings. Nb represent number of the rollers. ξ=3/2 for the rolling ball bearing used in the intermediate support. θj is the rotation angle of thejth roller at time t. γ and ωc are bearing radial clearance and rotational speed of the bearing retainer, respectively. uθj is the elastic radial deformation of the jth roller. r and R represent radius of the inner and outer ring. ωin and ωout 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(15)Mu¨+Gu˙+Ku=Fnonl+Funb.Funb is the unbalance force vector. Fnonl is the nonlinear forces exerted by the SFDs and rolling bearing in this work. M, K, and G 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 as(16)MIIMIJMJIMJJu¨Iu¨J+GIIGIJGJIGJJu˙Iu˙J+KIIKIJKJIKJJuIuJ=FIunb0+0FJnonl,where uI can be interpreted as linear DOFs, with only unbalance forces considered, while uJ 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 is(17a)uIuJ=ϕkϕc0IqkuJ=TqkuJ,where ϕk is the mass normalized normal mode matrix,ϕc is the mass normalized constrained mode, I is unity matrix, qk is the normal mode coordinates, and T is the transformation matrix. ϕc is given as below:(17b)ϕc=-KII-1KIJ.In (16), let uJ=0 and neglecting the gyroscopic effects gives(17c)MIIu¨I+KIIuI=0.Solving the eigenvalue problem corresponding to (17c) gives the mass normalized normal mode matrix ϕk and the modal frequency Ωk.Substitute (17a) into (16) and premultiply with TT to obtain the reduced equations of motion:(18a)IM-IJM-JIM-JJq¨ku¨J+G-IIG-IJG-JIG-JJq˙ku˙J+K-II00K-JJqkuJ=ϕkTFIunbϕcTFIunb+0FJnonl,where(18b)M-JJ=MJJ+MJIϕc+ϕcTMIIϕc+MIJ,M-IJ=M-JI=ϕkTMIIϕc+MIJ,G-II=ϕkTGIIϕk,G-JJ=GJJ+GJIϕc+ϕcTGIIϕc+GIJ+cJJ,G-IJ=G-JI=ϕkTGIIϕc+ϕkTGIJ,K-II=diag⁡Ωr2,1≤r≤n,K-JJ=KIJ+KJIϕc+kJJ.Ωr is the rth modal frequency obtained from (17c). FIunb is the unbalance force acting on linear DOFs. FJnonl is the nonlinear force acting on nonlinear DOFs.And(19)FJnonl=-kJJ+FJSFDFJB,where FJSFD=f1xsfdf1ysfd⋯fNxsfdfNysfd is the nonlinear force vector exerted by SFDs which can be calculated with (13a), (13b), (13c), (13d), and (13e). FJB=fxbfyb-fxb-fyb is the nonlinear force vector exerted by rolling bearing which can be calculated with (14a), (14b), (14c), (14d), and (14e). kJJ and cJJ are given below:(20)kJJ=kJB000,cJJ=cJB000,kJB=diag⁡kJiB,1≤i≤N,cJB=diag⁡cJiB,1≤i≤N,kJiB=kJi00kJi,cJiB=cJi00cJi.kJi and cJi are the stiffness and damping coefficients of the elastic support, respectively. N is the number of elastic supports.Rearranging (18a) and (18b) gives(21a)q¨k+G-IIq˙k+K-IIqk=ϕkTFIunb-M-IJu¨J-G-IJu˙J,(21b)M-JJu¨J+G-JJu˙J+K-JJuJ=FJnonl+ϕcTFIunb-M-JIq¨k-G-JIq˙k.In (21a) and (21b), the vectors qk and uJ 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 intervaltntn+1 is(22a)q˙kn+1u˙Jn+1=q˙knu˙Jn+1-βq¨knu¨Jn+βq¨kn+1u¨Jn+1Δt,(22b)qkn+1uJn+1=qknuJn+q˙knu˙JnΔt+0.5-αq¨knu¨Jn+αq¨kn+1u¨Jn+1Δt2,where tn+1=tn+Δt and Δt is the time increment. The superscripts n and n+1 denote tn and tn+1.The following equations can be obtained from (22a) and (22b):(23)q¨kn+1=aqkn+1-Aqn,q˙kn+1=bqkn+1-Bqn,(24)u¨Jn+1=auJn+1-AJn,u˙Jn+1=buJn+1-BJn,with(25)Aqn=1αΔt2qkn+1αΔtq˙kn+12α-1q¨kn,AJn=1αΔt2ukn+1αΔtu˙kn+12α-1u¨kn,Bqn=βαΔtqkn+βα-1q˙kn+βα-2q¨kn,BJn=βαΔtukn+βα-1u˙kn+βα-2u¨kn,a=1αΔt2,b=βαΔt.Substituting (22a) and (22b)–(25) into (21a) and (21b) yields(26a)qkn+1=Sq-1ϕkTFIunbn+1-VquJn+1+Wq,(26b)SJ-VJSq-1VquJn+1=FJnonln+1+ϕcT-VJSq-1ϕkTFIunbn+1-VJSq-1Wq+WJ,with(27a)Wq=M-IIAqn+G-IIBqn+M-IJAJn+G-IJBJn,(27b)WJ=M-JJAJn+G-JJBJn+M-JIAqn+G-JIBqn,(27c)Sq=aM-II+bG-II+K-II,(27d)SJ=aM-JJ+bG-JJ+K-JJ,(27e)Vq=aM-IJ+bG-IJ,(27f)VJ=aM-JI+bG-JI,(27g)FJnonln+1=FJnonlnuJn+1,u˙Jn+1.Use (24) to substitute for u˙Jn+1 in (27g) and substitute the new (27g) into (26b) to yield nonlinear equations of uJn+1. After solving (26b) for uJn+1 with numerical algorithms, qkn+1 can be obtained with uJn+1 substituted into (26a). Thus u˙Jn+1, u¨Jn+1, q˙kn+1, and q¨kn+1 can be solved from (23) and (24). The process continues to repeat and move forward to find uJn+2 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 Figure2.Figure 2
Flowchart of modeling and solving.


## 2.1. Finite Element Formulation

Figure1 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 as(1)qe=ue1ve1θe1φe1ue2ve2θe2φe2T.The displacement vector within the element can be interpolated as(2)ueξξ,tveξξ,tθeξξ,tφeξξ,t=N100N2N300N40N1-N200N3-N400-Nˇ1Nˇ200-Nˇ3Nˇ40Nˇ100Nˇ2Nˇ300Nˇ4qe.The shape functions in (2) are(3a)N1=1-1le2+12g12gleξ+3ξ2-2leξ3,N2=1le2+12gle2+6gξ-2le2+6gleξ2+ξ3,N3=1le2+12g12gleξ+3ξ2-2leξ3,N4=1le2+12g-6gξ-le2-6gleξ2+ξ3,(3b)Nˇ1=1le2+12g6leξ2-6ξ,Nˇ2=1le2+12gle2+12g-4le2+12gleξ+3ξ2,Nˇ3=1le2+12g6ξ-6leξ2,Nˇ4=1le2+12g12g-2le2leξ+3ξ2.In (3a) and (3b), le is the element length and g is formulated as given below:(4)g=EeIeκGeAe,where Ee is Young’s Modulus, Ie is the area moment of inertia, Ge is the shear modulus, and Ae is the cross-sectional area. κ denotes the shear coefficient, and for a hollow circular shaft section [16, 17](5)κ=61+μ1+λ227+12μ+4μ21+λ22+45+6μ+2μ2λ2,where μ is Poisson’s ratio and λ is the ratio of the inner shaft radius to the outer shaft radius. Hence for a solid shaft λ→0.Figure 1
Diagram of beam element.With rotary inertia and shear effects considered, the kinetic energy and strain energy for a shaft element are(6a)Te=12∫0leρeAeu˙eξ2ξ,t+ρeIeφ˙eξ2ξ,tdξ,(6b)Ue=12∫0le∂φeξξ,t∂ξTEeIe∂φeξξ,t∂ξdξ+12κ∫0leγTGeAeγdξ.And thus the element mass matrix is(7)me=ρeAe∫0leNTNdξ+ρeIe∫0legN′′′+N′TgN′′′+N′dξ.The element stiffness matrix is(8)ke=EeIe∫0leN′′TN′′dξ+EeIeg∫0leN′′′TN′′′dξ,where(9)N′=dN1N2N3N4dξ,N′′=d2N1N2N3N4dξ2,N′′′=d3N1N2N3N4dξ3.The element gyroscopic matrix is(10)ge=-2ρeIe∫0leB2TB1-B1TB2dξ,where(11)B1B2=0-Nˇ1Nˇ200-Nˇ3Nˇ40Nˇ100Nˇ2Nˇ300Nˇ4.For disk element, the mass and gyroscopic matrices are(12a)med=mdi0000mdi0000Idi0000Idi,(12b)ged=00000000000Ipi00-Ipi0,where mdi, Idi, and Ipi 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](13a)fxsfd=-μsRL3c2x2+y21/2xε˙I2+εψ˙I1-yε˙I1+εψ˙I3,fysfd=-μsRL3c2x2+y21/2yε˙I2+εψ˙I1+xε˙I1+εψ˙I3,(13b)ε=x2+y2c,(13c)ε˙=xx˙+yy˙cx2+y2,(13d)ψ˙=yx˙-xy˙x2+y2,(13e)tan⁡ψ=xy.x and y are the horizontal and vertical displacements of the journal. Ijj=1,2,3 are Sommerfeld integrals. R and L are radius and length of the SFD, respectively. μs and c 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] is(14a)fxb=kn∑j=1NbuθjξHuθjsin⁡θj,fyb=kn∑j=1NbuθjξHuθjcos⁡θj,(14b)θj=2πj-1Nb+ωct,(14c)uθj=xir-xorcos⁡θj+yir-yorsin⁡θj-γ2,(14d)Huθj=0,uθj≤0,uθj,uθj>0,(14e)ωc=ωinr+ωoutRR+r,where the superscripts ir and or in (14c) denote inner ring and outer ring of the bearing. kn in (14a) is the contact stiffness between rollers and rings. Nb represent number of the rollers. ξ=3/2 for the rolling ball bearing used in the intermediate support. θj is the rotation angle of thejth roller at time t. γ and ωc are bearing radial clearance and rotational speed of the bearing retainer, respectively. uθj is the elastic radial deformation of the jth roller. r and R represent radius of the inner and outer ring. ωin and ωout 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(15)Mu¨+Gu˙+Ku=Fnonl+Funb.Funb is the unbalance force vector. Fnonl is the nonlinear forces exerted by the SFDs and rolling bearing in this work. M, K, and G 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 as(16)MIIMIJMJIMJJu¨Iu¨J+GIIGIJGJIGJJu˙Iu˙J+KIIKIJKJIKJJuIuJ=FIunb0+0FJnonl,where uI can be interpreted as linear DOFs, with only unbalance forces considered, while uJ 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 is(17a)uIuJ=ϕkϕc0IqkuJ=TqkuJ,where ϕk is the mass normalized normal mode matrix,ϕc is the mass normalized constrained mode, I is unity matrix, qk is the normal mode coordinates, and T is the transformation matrix. ϕc is given as below:(17b)ϕc=-KII-1KIJ.In (16), let uJ=0 and neglecting the gyroscopic effects gives(17c)MIIu¨I+KIIuI=0.Solving the eigenvalue problem corresponding to (17c) gives the mass normalized normal mode matrix ϕk and the modal frequency Ωk.Substitute (17a) into (16) and premultiply with TT to obtain the reduced equations of motion:(18a)IM-IJM-JIM-JJq¨ku¨J+G-IIG-IJG-JIG-JJq˙ku˙J+K-II00K-JJqkuJ=ϕkTFIunbϕcTFIunb+0FJnonl,where(18b)M-JJ=MJJ+MJIϕc+ϕcTMIIϕc+MIJ,M-IJ=M-JI=ϕkTMIIϕc+MIJ,G-II=ϕkTGIIϕk,G-JJ=GJJ+GJIϕc+ϕcTGIIϕc+GIJ+cJJ,G-IJ=G-JI=ϕkTGIIϕc+ϕkTGIJ,K-II=diag⁡Ωr2,1≤r≤n,K-JJ=KIJ+KJIϕc+kJJ.Ωr is the rth modal frequency obtained from (17c). FIunb is the unbalance force acting on linear DOFs. FJnonl is the nonlinear force acting on nonlinear DOFs.And(19)FJnonl=-kJJ+FJSFDFJB,where FJSFD=f1xsfdf1ysfd⋯fNxsfdfNysfd is the nonlinear force vector exerted by SFDs which can be calculated with (13a), (13b), (13c), (13d), and (13e). FJB=fxbfyb-fxb-fyb is the nonlinear force vector exerted by rolling bearing which can be calculated with (14a), (14b), (14c), (14d), and (14e). kJJ and cJJ are given below:(20)kJJ=kJB000,cJJ=cJB000,kJB=diag⁡kJiB,1≤i≤N,cJB=diag⁡cJiB,1≤i≤N,kJiB=kJi00kJi,cJiB=cJi00cJi.kJi and cJi are the stiffness and damping coefficients of the elastic support, respectively. N is the number of elastic supports.Rearranging (18a) and (18b) gives(21a)q¨k+G-IIq˙k+K-IIqk=ϕkTFIunb-M-IJu¨J-G-IJu˙J,(21b)M-JJu¨J+G-JJu˙J+K-JJuJ=FJnonl+ϕcTFIunb-M-JIq¨k-G-JIq˙k.In (21a) and (21b), the vectors qk and uJ 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 intervaltntn+1 is(22a)q˙kn+1u˙Jn+1=q˙knu˙Jn+1-βq¨knu¨Jn+βq¨kn+1u¨Jn+1Δt,(22b)qkn+1uJn+1=qknuJn+q˙knu˙JnΔt+0.5-αq¨knu¨Jn+αq¨kn+1u¨Jn+1Δt2,where tn+1=tn+Δt and Δt is the time increment. The superscripts n and n+1 denote tn and tn+1.The following equations can be obtained from (22a) and (22b):(23)q¨kn+1=aqkn+1-Aqn,q˙kn+1=bqkn+1-Bqn,(24)u¨Jn+1=auJn+1-AJn,u˙Jn+1=buJn+1-BJn,with(25)Aqn=1αΔt2qkn+1αΔtq˙kn+12α-1q¨kn,AJn=1αΔt2ukn+1αΔtu˙kn+12α-1u¨kn,Bqn=βαΔtqkn+βα-1q˙kn+βα-2q¨kn,BJn=βαΔtukn+βα-1u˙kn+βα-2u¨kn,a=1αΔt2,b=βαΔt.Substituting (22a) and (22b)–(25) into (21a) and (21b) yields(26a)qkn+1=Sq-1ϕkTFIunbn+1-VquJn+1+Wq,(26b)SJ-VJSq-1VquJn+1=FJnonln+1+ϕcT-VJSq-1ϕkTFIunbn+1-VJSq-1Wq+WJ,with(27a)Wq=M-IIAqn+G-IIBqn+M-IJAJn+G-IJBJn,(27b)WJ=M-JJAJn+G-JJBJn+M-JIAqn+G-JIBqn,(27c)Sq=aM-II+bG-II+K-II,(27d)SJ=aM-JJ+bG-JJ+K-JJ,(27e)Vq=aM-IJ+bG-IJ,(27f)VJ=aM-JI+bG-JI,(27g)FJnonln+1=FJnonlnuJn+1,u˙Jn+1.Use (24) to substitute for u˙Jn+1 in (27g) and substitute the new (27g) into (26b) to yield nonlinear equations of uJn+1. After solving (26b) for uJn+1 with numerical algorithms, qkn+1 can be obtained with uJn+1 substituted into (26a). Thus u˙Jn+1, u¨Jn+1, q˙kn+1, and q¨kn+1 can be solved from (23) and (24). The process continues to repeat and move forward to find uJn+2 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 Figure2.Figure 2
Flowchart of modeling and solving.


## 3. Test Rig Description

The coaxial test rig studied in this paper is presented in Figure3.Figure 3
Structural diagram of the coaxial rotor system.The stationary coordinate system in Figure3 consists of three mutually perpendicular axes, ox, oy, and oz, intersecting at the point o and axis of the rotor coincides with axis oz. The axes ox and oz are horizontal, and oy is vertical.The rotor system consists of two shafts disposed along the same axisoz, 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 Tables1–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 ϕk. For the Newmark-beta method, α=0.25 and β=0.5.Table 1
Dimension and elements information of the rotor system.
Node number
Axial location (m)
Bearing/disk
Element number
Outer diameter (m)
Inner diameter (m)
1
0
1
0.018
0.00
2
0.08143
2
0.018
0.00
3
0.16286
3
0.018
0.00
4
0.24429
4
0.018
0.00
5
0.24909
5
0.018
0.00
6
0.25479
6
0.018
0.00
7
0.28879
7
0.018
0.00
8
0.32279
8
0.018
0.00
9
0.35879
Disk number 1
9
0.018
0.00
10
0.38369
10
0.018
0.00
11
0.40859
11
0.018
0.00
12
0.43349
12
0.018
0.00
13
0.43869
Bearing number 1
13
0.018
0.00
14
0.44479
14
0.022
0.00
15
0.54752
15
0.022
0.00
16
0.65025
16
0.022
0.00
17
0.75298
17
0.022
0.00
18
0.85571
18
0.022
0.00
19
0.95844
19
0.022
0.00
20
1.06117
20
0.022
0.00
21
1.06517
Bearing number 4
21
0.022
0.00
22
1.06867
22
0.022
0.00
23
1.08867
23
0.022
0.00
24
1.10867
Disk number 2
24
0.022
0.00
25
1.14274
25
0.022
0.00
26
1.17681
26
0.022
0.00
27
1.21088
27
0.017
0.00
28
1.21488
Bearing number 2
28
0.014
0.00
29
1.21838
29
0.014
0.00
30
1.23038
End of inner rotor
31
0.64200
30
0.035
0.03
32
0.66065
31
0.035
0.03
33
0.67930
Bearing number 3
32
0.035
0.03
34
0.68650
33
0.038
0.03
35
0.71170
34
0.038
0.03
36
0.73690
35
0.038
0.03
37
0.76210
Disk number 3
36
0.038
0.03
38
0.80784
37
0.038
0.03
39
0.85358
38
0.038
0.03
40
0.89932
39
0.038
0.03
41
0.94506
40
0.038
0.03
42
0.99080
Disk number 4
41
0.038
0.03
43
1.01430
42
0.038
0.03
44
1.02030
43
0.070
0.03
45
1.03330
44
0.060
0.03
46
1.06030
45
0.060
0.03
47
1.06430
Bearing number 4
46
0.060
0.03
48
1.07380Table 2
Stiffness of elastic supports (squirrel cage).
Support I
Support II
Support III
Support IV
Stiffness (N/m)
1.45 × 106
2.21 × 105
9.29 × 105
—Table 3
Parameters of the intermediate bearing.
Radius of inner ring (mm)
Radius of outer ring (mm)
Number of rollers
Contact stiffness (N/m3/2)
Radial clearance (μm)
9.37
14.13
9
7.055 × 109
6Table 4
Parameters of SFDs.
Inner rotor
Outer rotor
SFD I
SFD II
SFD III
Radius R/mm
25
18
35
Length L/mm
15
15
20
Radial clearance c/mm
0.1
0.1
0.08
Dynamic viscosityμs/10−2 Pa⋅s
1.0752Table 5
Unbalance configuration and inertia properties of disks.
Inner rotor
Outer rotor
Disk 1
Disk 2
Disk 3
Disk 4
Unbalance (×10−5 kg⋅m)
2
4
1
2
Mass (kg)
2.3386
2.3386
3.2590
1.6303
Polar moment of inertia (kg⋅m2)
0.00815
0.00815
0.01561
0.00661Geometric dimensions and information of each element are listed in Table1. 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λ=ω2/ω1, where ω2 is the rotational speed of the outer rotor and ω1 is the rotational speed of the inner rotor. λ=1.65 (corotating) and λ=-1.65 (counterrotating) cases are studied in the current work.


### 4.1. Overview of the Results

With the model established by method described in Section2, 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 Figures4 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,ω1 and ω2, 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ω1 and ω2 in the responses of disks 2 and 4, such as 2ω1+ω2, ω1+2ω2, 3ω1+2ω2, 2ω1+3ω2, and (ω2-ω1)/2, 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 2ω1+ω2, ω1+2ω2, 3ω1+2ω2, and 2ω1+3ω2, while in Figure 5(a) only 2ω1+ω2 and ω1+2ω2 are observed. In Figure 4(b), the combination frequency components are 2ω1+ω2, ω1+2ω2, 3ω1+2ω2, 2ω1+3ω2, and (ω2-ω1)/2 while in Figure 5(b) those combinations are 2ω2-ω1,3ω2-2ω1, and6ω2+ω1/3.(3)
In Figure4, 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.Figure 4
Spectrum cascade of the horizontal response of disks 2 and 4, numerical simulation.λ=-1.65.
(a)
Disk 2
(b)
Disk 4Figure 5
Spectrum cascade of the horizontal response of disks 2 and 4, numerical simulation.λ=1.65.
(a)
Disk 2
(b)
Disk 4Figure 6
Spectrum cascade of the horizontal response of disks 2 and 4, experimental results.λ=-1.65.
(a)
Disk 2
(b)
Disk 4Figure 7
Spectrum cascade of the horizontal response of disks 2 and 4, experimental results.λ=1.65.
(a)
Disk 2
(b)
Disk 4


### 4.2. Unbalance Response Analysis

Respectively, for counter- and corotating cases, Figures8 and 9 show bifurcation diagrams for horizontal displacement of disk 2 and disk 4 with rotor speed. Sampling period for the bifurcation diagram is 5×2π/ω1 because the rotational speed ratio of outer and inner rotor is ±1.65 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.Figure 8
Bifurcation diagrams of the horizontal response of disk 2 and disk 4 with rotor speed.λ=-1.65.
(a)
Disk 2
(b)
Disk 4
(c)
Disk 2
(d)
Disk 4Figure 9
Bifurcation diagrams of the horizontal response of disk 2 and disk 4 with rotor speed.λ=1.65.
(a)
Disk 2
(b)
Disk 4
(c)
Disk 2
(d)
Disk 4Basically, disks 2 and 4 execute multiple period orbital motion around critical speeds, which can be seen from the bifurcation diagrams shown in Figures8 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 ω1 and ω2 are observed in the spectrogram shown in Figures 10(a) and 11(a). And the amplitude of ω2 is much larger than that of ω1, which means the outer rotor is the main source of vibration.Figure 10
Response analysis of disk 2:ω1=102rad/s and λ=-1.65.
(a)
Spectrogram
(b)
Orbit
(c)
Phase trajectory
(d)
Poincare mapFigure 11
Response analysis of disk 2:ω1=102rad/s and λ=1.65.
(a)
Spectrogram
(b)
Orbit
(c)
Phase trajectory
(d)
Poincare mapAs 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, Figures12 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 ω1 and ω2 can be found for the conditions when no obvious periodicity can be observed.Figure 12
Response analysis of disk 2:ω1=114rad/s and λ=-1.65.
(a)
Spectrogram
(b)
Orbit
(c)
Phase trajectory
(d)
Poincare mapFigure 13
Response analysis of disk 2:ω1=114rad/s and λ=1.65.
(a)
Spectrogram
(b)
Orbit
(c)
Phase trajectory
(d)
Poincare mapFigure 14
Response analysis of disk 2:ω1=140rad/s and λ=-1.65.
(a)
Spectrogram
(b)
Orbit
(c)
Phase trajectory
(d)
Poincare mapFigure 15
Response analysis of disk 2:ω1=140rad/s and λ=1.65.
(a)
Spectrogram
(b)
Orbit
(c)
Phase trajectory
(d)
Poincare mapFigure 16
Response analysis of disk 2:ω1=170rad/s and λ=-1.65.
(a)
Spectrogram
(b)
Orbit
(c)
Phase trajectory
(d)
Poincare mapFigure 17
Response analysis of disk 2:ω1=170rad/s and λ=1.65.
(a)
Spectrogram
(b)
Orbit
(c)
Phase trajectory
(d)
Poincare mapAt 170 rad/s, approximately the next critical speed, four-period quasi-periodic motion emerges again; see Figures16 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, Figures10, 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, Figures4–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.Figure 18
Orbit of disk 2 and disk 4 withω1=128rad/s and λ=-1.65.
(a)
Disk 2—numerical results
(b)
Disk 2—experimental results
(c)
Disk 4—numerical results
(d)
Disk 4—experimental resultsFigure 19
Orbit of disk 2 and disk 4 withω1=160rad/s and λ=-1.65.
(a)
Disk 2—numerical results
(b)
Disk 2—experimental results
(c)
Disk 4—numerical results
(d)
Disk 4—experimental resultsFigure 20
Orbit of disk 2 and disk 4 withω1=210rad/s and λ=-1.65.
(a)
Disk 2—numerical results
(b)
Disk 2—experimental results
(c)
Disk 4—numerical results
(d)
Disk 4—experimental resultsFigure 21
Orbit of disk 2 and disk 4 withω1=146rad/s and λ=1.65.
(a)
Disk 2—numerical results
(b)
Disk 2—experimental results
(c)
Disk 4—numerical results
(d)
Disk 4—experimental resultsFigure 22
Orbit of disk 2 and disk 4 withω1=170rad/s and λ=1.65.
(a)
Disk 2—numerical results
(b)
Disk 2—experimental results
(c)
Disk 4—numerical results
(d)
Disk 4—experimental results


## 4.1. Overview of the Results

With the model established by method described in Section2, 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 Figures4 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,ω1 and ω2, 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ω1 and ω2 in the responses of disks 2 and 4, such as 2ω1+ω2, ω1+2ω2, 3ω1+2ω2, 2ω1+3ω2, and (ω2-ω1)/2, 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 2ω1+ω2, ω1+2ω2, 3ω1+2ω2, and 2ω1+3ω2, while in Figure 5(a) only 2ω1+ω2 and ω1+2ω2 are observed. In Figure 4(b), the combination frequency components are 2ω1+ω2, ω1+2ω2, 3ω1+2ω2, 2ω1+3ω2, and (ω2-ω1)/2 while in Figure 5(b) those combinations are 2ω2-ω1,3ω2-2ω1, and6ω2+ω1/3.(3)
In Figure4, 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.Figure 4
Spectrum cascade of the horizontal response of disks 2 and 4, numerical simulation.λ=-1.65.
(a)
Disk 2
(b)
Disk 4Figure 5
Spectrum cascade of the horizontal response of disks 2 and 4, numerical simulation.λ=1.65.
(a)
Disk 2
(b)
Disk 4Figure 6
Spectrum cascade of the horizontal response of disks 2 and 4, experimental results.λ=-1.65.
(a)
Disk 2
(b)
Disk 4Figure 7
Spectrum cascade of the horizontal response of disks 2 and 4, experimental results.λ=1.65.
(a)
Disk 2
(b)
Disk 4


## 4.2. Unbalance Response Analysis

Respectively, for counter- and corotating cases, Figures8 and 9 show bifurcation diagrams for horizontal displacement of disk 2 and disk 4 with rotor speed. Sampling period for the bifurcation diagram is 5×2π/ω1 because the rotational speed ratio of outer and inner rotor is ±1.65 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.Figure 8
Bifurcation diagrams of the horizontal response of disk 2 and disk 4 with rotor speed.λ=-1.65.
(a)
Disk 2
(b)
Disk 4
(c)
Disk 2
(d)
Disk 4Figure 9
Bifurcation diagrams of the horizontal response of disk 2 and disk 4 with rotor speed.λ=1.65.
(a)
Disk 2
(b)
Disk 4
(c)
Disk 2
(d)
Disk 4Basically, disks 2 and 4 execute multiple period orbital motion around critical speeds, which can be seen from the bifurcation diagrams shown in Figures8 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 ω1 and ω2 are observed in the spectrogram shown in Figures 10(a) and 11(a). And the amplitude of ω2 is much larger than that of ω1, which means the outer rotor is the main source of vibration.Figure 10
Response analysis of disk 2:ω1=102rad/s and λ=-1.65.
(a)
Spectrogram
(b)
Orbit
(c)
Phase trajectory
(d)
Poincare mapFigure 11
Response analysis of disk 2:ω1=102rad/s and λ=1.65.
(a)
Spectrogram
(b)
Orbit
(c)
Phase trajectory
(d)
Poincare mapAs 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, Figures12 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 ω1 and ω2 can be found for the conditions when no obvious periodicity can be observed.Figure 12
Response analysis of disk 2:ω1=114rad/s and λ=-1.65.
(a)
Spectrogram
(b)
Orbit
(c)
Phase trajectory
(d)
Poincare mapFigure 13
Response analysis of disk 2:ω1=114rad/s and λ=1.65.
(a)
Spectrogram
(b)
Orbit
(c)
Phase trajectory
(d)
Poincare mapFigure 14
Response analysis of disk 2:ω1=140rad/s and λ=-1.65.
(a)
Spectrogram
(b)
Orbit
(c)
Phase trajectory
(d)
Poincare mapFigure 15
Response analysis of disk 2:ω1=140rad/s and λ=1.65.
(a)
Spectrogram
(b)
Orbit
(c)
Phase trajectory
(d)
Poincare mapFigure 16
Response analysis of disk 2:ω1=170rad/s and λ=-1.65.
(a)
Spectrogram
(b)
Orbit
(c)
Phase trajectory
(d)
Poincare mapFigure 17
Response analysis of disk 2:ω1=170rad/s and λ=1.65.
(a)
Spectrogram
(b)
Orbit
(c)
Phase trajectory
(d)
Poincare mapAt 170 rad/s, approximately the next critical speed, four-period quasi-periodic motion emerges again; see Figures16 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, Figures10, 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, Figures4–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.Figure 18
Orbit of disk 2 and disk 4 withω1=128rad/s and λ=-1.65.
(a)
Disk 2—numerical results
(b)
Disk 2—experimental results
(c)
Disk 4—numerical results
(d)
Disk 4—experimental resultsFigure 19
Orbit of disk 2 and disk 4 withω1=160rad/s and λ=-1.65.
(a)
Disk 2—numerical results
(b)
Disk 2—experimental results
(c)
Disk 4—numerical results
(d)
Disk 4—experimental resultsFigure 20
Orbit of disk 2 and disk 4 withω1=210rad/s and λ=-1.65.
(a)
Disk 2—numerical results
(b)
Disk 2—experimental results
(c)
Disk 4—numerical results
(d)
Disk 4—experimental resultsFigure 21
Orbit of disk 2 and disk 4 withω1=146rad/s and λ=1.65.
(a)
Disk 2—numerical results
(b)
Disk 2—experimental results
(c)
Disk 4—numerical results
(d)
Disk 4—experimental resultsFigure 22
Orbit of disk 2 and disk 4 withω1=170rad/s and λ=1.65.
(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.



---

*Source: 5493763-2017-06-07.xml*