Dynamic modeling and analysis of the rotor–stator coupling system of a coaxial contra-rotating gearbox

In this paper, a main gearbox using an encased differential gear train to achieve coaxial contra-rotating is considered, and a dynamic modeling method of the rotor–stator coupling system of the gearbox based on box model updating is introduced. The transverse torsional dynamic model of the gear transmission subsystem is established based on the lumped parameter method. The finite element model of the box is updated according to the modal test data, and the reduced dynamic parameters of the box are obtained. According to the displacement coordination condition, the dynamic model of the rotor–stator system of the gearbox is established. The vibration response of the transmission components with or without the coupling box is calculated by numerical integration, and the response of the box caused by the dynamic support reaction force is analyzed by the finite element method. The results show that the vibration peak and fluctuation range of the transmission parts with coupling box are smaller than those without coupling box. The box response at the support of the input bevel gear pair is large, while that at the support of the output shaft is small.


Dynamic model of the transmission system
The transmission system of the coaxial main gearbox consists of a spiral bevel gear pair and an encased differential gear train, as shown in Fig. 1.The encased differential gear train can be divided into two parts: the encased stage composed of sun gear s 1 , stepped planets a i , b i and ring gear r 1 , and the differential stage composed of sun gear s 2 , planet p j , ring gear r 2 and carrier c 2 .
When the input power is transmitted to the encased differential gear train through the spiral bevel gear pair, it is divided into two paths: one part is transmitted from sun gear s 1 to ring gear r 1 , and the other part is transmitted from sun gear s 2 to carrier c 2 and ring gear r 2 through the differential stage.The transmission power of the carrier c 2 is output by the inner shaft c s , and the power of the ring gear r 1 and r 2 is output by the outer shaft r s after confluence.
Through the proper parameter design, as shown in Table 1, the inner and outer output shafts can rotate in opposite directions at the same speed.In the system, the number of planets p is N, and the number of stepped planets ab is M. The lumped parameter dynamic model of the differential stage gear train is established based on the fixed coordinate system OXY, whose origin is the theoretical installation center of sun gear s 2 , as shown in Fig. 2. The local coordinate system O pj ζ pj η pj is fixed on the carrier and rotates with it at a constant speed, and its origin is located at the theoretical installation center O pj of the j-th planet.
When the position coordinate of the centroid of planet p j in the local coordinate system is ( ζ pj , η pj ), the dis- tance OO ′ pj in the coordinate system OXY can be expressed as: According to the coordinate transformation rules, the absolute acceleration components of the planet p j in the O pj ζ pj and O pj η pj directions are shown as follows 21 : There are Without considering nonlinearities such as clearance, the dynamic model of the transmission system of the gearbox can be obtained by using the lumped parameter method, as shown in Fig. 3.The origin of the overall coordinate system of the transmission system is located at the intersection of the axes of the bevel gear pair.The positive direction of the X-axis of the coordinate system points to the input shaft, and the positive direction of the Z-axis points to the output shaft.
The dynamic model assumes that the size, mass, and moment of inertia of all planets in each planetary set are equal, the bearings of each planet have the same support stiffness and the stiffness in all directions is also equal.
The parameters in the model are defined as follows: k m is the mesh stiffness of gear pair m; k h (h = s 1 , r s , c s , a, b, p) is the radial support stiffness of member h; k h1x (k h2x ), k h1y (k h2y ) and k h1z (k h2z ) are the three-way support stiffness of bevel gear h 1 (h 2 ); and k b AB and k u AB are the radial coupling stiffness and torsional stiffness between members A and B, respectively.The definition of the damping symbol c is similar.
The vibration displacement vector of the system can be written as: The vector includes L = 29 + 3(2M + N) degrees of freedom, where x, y, and z are the vibration displacements in the translational direction of the bevel gear, center gear, carrier and output shaft.ζ and η are the vibration displacements of the planet along the radial and tangential directions of the carrier, respectively, and u is the torsional vibration line displacement of the component.
In the meshing cycle of a gear pair, the difference in mesh stiffness between single-tooth meshing and doubletooth meshing will cause large vibration acceleration of the gear system, which is an important excitation source of vibration of the gear system.The change in the mesh stiffness in the single-tooth meshing area or the doubletooth meshing area is small, which has little effect on the vibration response of the gear system 22 .Therefore, the time-varying mesh stiffness of gear pair m considering the mesh phase can be expanded into a Fourier series expression 23 : where where k a m and 2k v m are the mean value and fluctuation value of the mesh stiffness of gear pair m, respectively, and γ m is the mesh phase difference.
The mesh damping of gear pair m is: q= x s1 , y s1 , u s1 , x r1 , y r1 , u r1 , ζ ai , η ai , u ai , ζ bi , η bi , u bi , x s2 , y s2 , u s2 , x r2 ,y r2 , u r2 , x c2 , y c2 , u c2 , ζ pj , η pj , u pj , x cs , y cs , u cs , x rs , y rs , u rs ,x h1 , y h1 , z h1 , u h1 , x h2 , y h2 , z h2 , u h2 where ζ m is the mesh damping ratio of gear pair m, generally taken as 0.03-0.17;I m1 and I m2 are the moments of inertia of the driving and driven gears, respectively; and r bm1 and r bm2 are the base circle radii of the driving and driven gears, respectively.

Relative displacement
Compression is defined as positive, and the relative displacement of the sun gear and planet along the meshing line is: Similarly, the relative displacement between the ring gear and the planet is defined as: where α ′ r1 and α ′ s1 are the working pressure angles of the internal and external meshing pair of the encased stage, respectively; α ′ r2 and α ′ s2 are the angles of the differential stage; and ψ ai , ψ bi and ψ pj are the position angles of planet pj, which can be expressed as: The relative displacements of planet p j and carrier c 2 in the radial and tangential directions are as follows: The relative displacement of the bevel gear pair is: where In the formula, α h2 , δ h2 , and β h2 are the normal pressure angle, pitch cone angle, and helix angle at the mid- point of the tooth width of bevel gear h 2 , respectively.e m is the comprehensive error of gear pair m along the meshing line.

Dynamic differential equations
Based on the translational-torsional dynamics model, the differential equations of the transmission system are derived as follows: 1. Dynamic equations of encased stage sun gear s 1 (5)    (32)  www.nature.com/scientificreports/these frequencies have no significant regularity and are all within 5%.The reason for the error is the difference in geometric dimensions and material properties between the machined box and the finite element model.

Model updating based on test data
Based on the modal test data, the finite element model of the box can be updated, and the accurate dynamic parameters of the box can be extracted to establish the rotor-stator coupling dynamic model.Considering the manufacturing process of the box, the main geometric and material parameters were selected as variable parameters: diameter, height, wall thickness, elastic modulus, and Poisson's ratio, and the parametric finite element calculation model for the box is established.
The minimum error between the calculation and experimental modal data is defined as the objective function.Considering that the effective mass of low-order frequencies accounts for the vast majority of the actual mass, only the first six nonzero natural frequencies are considered in the update process.Then,   In the formula, f ti is the frequency measured by the modal test, f fi is the frequency calculated by the finite element, and v max j and v min j are the upper and lower limits of the variable parameter v j .The model updating is transformed into a constrained structural optimization problem, and an optimization model based on the response surface is established based on ANSYS software, as shown in Fig. 11.The error between the simulation frequency and the test frequency is defined as the output parameter, and the goal of minimizing the output parameter can be achieved by adjusting the variable parameters.

Updated results.
Correlation analysis of the parameters is carried out, and the response surface is calculated.After approximately 140 iterations, the target value remained stable and the updated parameters are obtained.Among them, the elastic modulus is 195 GPa, Poisson's ratio is 0.255, and density is 7830 kg/m 3 .Based on the updated material parameters and geometric parameters, the finite element model of the box is rebuilt, and the first six natural frequencies are obtained by modal analysis, as shown in Fig. 12.
The frequencies calculated by the updated finite element model and those obtained by the test are shown in Table 4.It can be seen from the table that the error of the first five order frequencies after model updating is significantly reduced, while the error of the 6th-order frequency increases slightly.

Box substructure modeling
Based on the updated finite element model, the center of the support bearing hole on the box is established as the main node, as shown in Fig. 13.The modal analysis of the box is carried out to dynamically condense the internal degrees of freedom of the box.The stiffness, mass, and damping matrix of the box are extracted at the main nodes to realize order reduction and condensation and finally reduce the degree of freedom.

Dynamic modeling and analysis of the rotor-stator coupling system of the gearbox
The center of the supporting bearing is defined as the coupling point, and the dynamic equation of the rotor-stator coupling system of the gearbox can be obtained by combining the two substructure equations according to the displacement coordination.where M t , C t , and K t are the mass, damping, and stiffness matrices of the transmission substructure, respectively.M c , C c , and K c are the mass, damping, and stiffness matrices of the box substructure.K tc (K ct ) and C tc (C ct ) are the coupling stiffness and damping matrix of the transmission substructure and box substructure, respectively.q t and q c are the vibration displacement vectors of the transmission substructure and box substructure, respectively.Q t and Q c are the load vectors of the transmission substructure and box substructure, respectively.Based on the equation of the coupling system, the vibration response characteristics of the transmission system considering the flexibility of the box can be analysed.

Vibration response of the transmission system
The load torque of the inner and outer shafts is taken as T r = T c = 2815 N m, and the input speed is changed.The vibration response of the transmission components is obtained by solving Eqs. ( 53) and ( 55) to analyze the influence of the flexibility of the box on the vibration characteristics of the transmission system.
The vibration displacement response of the input bevel gear with the rotation speed is shown in Fig. 14.It can be seen from the figure that with the increase in the input rotation speed, the vibration displacement of the input bevel gear in the system with or without coupling box has varying degrees of resonance peaks.
The vibration displacement of the rotor-stator coupling system is larger than that without the coupling box, and the fluctuation range is smaller than that without the coupling box.This is because the overall stiffness of the system is reduced after considering the flexibility of the box, and part of the vibration energy is absorbed, which makes the fluctuations more stable.
The vibration displacement of the outer output shaft is shown in Fig. 15.Similarly, the vibration displacement fluctuation ranges of the rotor-stator coupling system coupling the box are small, and the peak value is offset.the support of the input bevel gear pair and output shaft, as well as the installation position of the gearbox, are analyzed, as shown in Fig. 18.
The root mean square value of normal vibration acceleration at each response point obtained by simulation is shown in Table 5.
It can be seen from Figs. 1, 18 and Table 5 that the vibration response of the box near the support of the input bevel gear pair is the largest, which is due to the large excitation of the box caused by the support reaction force caused by the meshing force of the input stage bevel gear pair, while the vibration response at the support of the output shaft and the installation of the gearbox is small.

Conclusions
Based on the modal test data, the finite element model of the box is updated.The stiffness matrix, mass matrix, and damping matrix at the main node of the updated box are extracted using the super element method, achieving the reduction and condensation of the box model and reducing the degree of freedom.A coupled dynamic model of the coaxial contra-rotating gearbox rotor-stator system is established using displacement coordination conditions, and the vibration characteristics of the system are analyzed based on this model.The results of this study can be summarized as follows: 1.The optimization method based on the response surface can effectively update the finite element calculation model of the box and reduce the error between the finite element calculation and modal test.2. With increasing input speed, the vibration response of various components in the system will exhibit resonance peaks to varying degrees whether the box is considered or not.When considering the flexibility of the box, the fluctuation of the component is smaller than that without the box.3. The results of the vibration response of the box calculated by the transient structure method show that the response of the support of the bevel gear pair is large and that of the output shaft is small.

Figure 1 .
Figure 1.Transmission system of coaxial contra-rotating main gearbox.

Figure 3 .
Figure 3. Translational-torsional dynamic model of the transmission system.

Figure 4 .
Figure 4. Finite element model of the box.

Figure 5 .
Figure 5. Natural frequencies of the box.

Figure 6 .
Figure 6.Vibration modes of the box.

Figure 9 .
Figure 9. Natural frequencies of the box.

Figure 10 .
Figure 10.Vibration mode of the box.

Figure 11 .
Figure 11.Optimization model based on the response surface.

Figure 12 .
Figure 12.Natural frequencies of the updated box.

Figure 13 .
Figure 13.Finite element polycondensation model of the box.

Figure 14 .
Figure 14.Vibration displacement of the input bevel gear.

Figure 15 .
Figure 15.Vibration displacement of the outer output shaft.

Figure 18 .
Figure 18.Vibration response analysis model of the box.

Table 1 .
The parameters of the gearbox.

Table 3 .
Frequencies of modal calculation and testing.

Table 4 .
Frequencies of the updated modal calculation and modal test.