An Improved Computational Method for Vibration Response and Radiation Noise Analysis of Two-Stage Gearbox

In this work, an improved computational method for the prediction of the vibration and noise of a gearbox that considers the flexibility of the shaft is developed. Based on the finite element method (FEM), a coupled dynamic model of a spur gear-shaft-bearing system is established, and the time-varying mesh stiffness (TVMS), the time-varying bearing stiffness (TVBS), and the flexibility of the shaft are considered. The Newmark integration method (NIM) is utilized to obtain the dynamic load of the bearing. Furthermore, the proposed model is validated by experiments. The bearing load is then considered to be the excitation of the housing, and the radiated noise is calculated via the finite element method/boundary element method (FEM/BEM). The effects of the shaft flexibility on the bearing response and radiated noise are discussed based on the proposed method. The results demonstrate that, when the shaft flexibility is considered, the system undergoes the bending vibration of the shaft, and the vibration amplitude and excitation frequency components of the bearing load decrease significantly. Additionally, the main resonance mode of the gearbox is changed, and the radiated noise is enhanced. The effects of the input speed and shaft stiffness on the bearing response and radiated noise are also investigated. The results provide a theoretical basis for the further development of the vibration and noise reduction of gearboxes.


I. INTRODUCTION
Gear systems are characterized by compact structure and high transmission efficiency, and they are widely used in the automotive, aerospace and power generation industries. Higher reliability and lower noise are desirable for modern gear transmissions. Gearbox vibration originates mostly from the meshing action of the gear teeth, and is then transmitted through the shafts and bearings to the housing and surrounding environment. The strong vibration of the gearbox is not only adverse to stability but also radiates noise to the outer space, resulting in noise pollution and affecting the equipment performance. For the reason, it is of great significance to estimate and control the vibration and noise of gearboxes to improve the dynamic characteristics of gear systems and design a high-quality silent gearbox.
Previous studies have yielded a remarkable variety of models used in gear system dynamics. The lumped-parameters The associate editor coordinating the review of this manuscript and approving it for publication was Hassen Ouakad . model is widely used in the numerical simulation of the dynamic behaviors of gear systems [1]- [3]; in this type of model, the gear, shaft, and bearing are considered as one, and are represented by mass blocks, springs, and damping elements. However, the influences of the shaft structure and mass distribution are neglected; thus, this model cannot accurately predict the dynamic responses of the shafts and bearings. In contrast, the finite element model has good calculation accuracy and is of greater importance for the prediction of the vibration and noise of gearboxes. Neriva et al. [4] applied finite element theory to study the gear rotor system, and considered the elastic shafts to be equivalent to the discrete stiffness and mass matrix. On this basis, Kahraman et al. [5] established a dynamic model of a gearshaft-bearing system to study the responses of gears and bearings. Later, Kubur et al. [6] proposed a dynamic model of a multi-shaft helical gear reduction unit formed by N flexible shafts, and demonstrated the influences of a number of system parameters including shaft angles, shaft dimensions, and bearing stiffness. Subsequently, Chang et al. [7] VOLUME 8,2020 This work is licensed under a Creative Commons Attribution 4.0 License. For more information, see https://creativecommons.org/licenses/by/4.0/ proposed a comprehensive, fully-coupled dynamic model of a gear-shaft-bearing-case system to determine more accurate bearing responses to predict the noise of the gearbox. Saxena et al. [8] studied the effects of the mesh stiffness and damping due to gear pairs on the modal characteristics of a geared flexible rotor-shaft system using Timoshenko beam elements. Xu et al. [9] further proposed an improved dynamic model of gearbox that considers flexibility of shaft and housing and suitable for variable speed process analysis. Yuan et al. [10] developed a quasi-static contact analysis model for a wide-faced cylindrical geared rotor system that considers shaft deflection, and investigated the effects of the supporting layout, power transmission path and output torque on the quasi-static and dynamic behaviors of the system. In the early stages of research, empirical formulas for the estimation of the noise intensity of a gear pair were obtained based on many experiments [11]- [13]. However, these acoustic estimates may not be suitable for other gearboxes, as the radiated noise depends strongly on the shape of the housing. The finite element method/boundary element method (FEM/BEM) has since been widely applied to the numerical analysis of the reducer. Kato et al. [14] simulated the vibration and noise radiation of a single-stage gearbox by combining the finite element analysis of the vibration and the boundary element analysis of the noise. The results of this numerical analysis were compared with experimental results. Later, Kostic and Ognjanovic [15] found that the noise emission of gearboxes depends on the disturbances, insulating capabilities and modal behavior of the housing. In recent years, investigations of the factors that affect these radiated noises have started to constitute a prominent area of research in the field. Abbes et al. [16] developed a finite element formulation for the modeling of the dynamic behavior of an interior coupled acoustic-structural system and studied the effect of the fluid inside the gearbox housing. Subsequently, Zhou et al. [17] and Zhou [18] presented a linear vibration model using the lumped-mass method, and FEM/BEM was used to analyze the vibration and noise of the reducer under various operation conditions. Guo et al. [19] developed a system-level vibro-acoustic model of an actual gearbox and compared the numerical results to experimental results. The influences of a standard rolling element and modified journal bearings on the radiated noise of a gearbox were investigated. Xi et al. [20] investigated the effects of different macro parameters, including the pressure angle and helical angle, on the vibration and noise of a gear transmission system under the rated working condition.
Furthermore, numerous scholars have focused on optimizing the radiated noise of the gearbox in different ways. Tanaka et al. [21] developed a method for obtaining the distribution of the sound pressure and identifying the areas from which noise radiates intensely, and designed a low-noise gearbox via the addition of ribs. Ene and Dimofte [22] measured the gear mesh noise with a microphone suspended above the gearbox, and found that wave bearings can attenuate the gear mesh noise by more than 10 dB as compared to rolling element bearings. Jolivet et al. [23] found that finishing processes can be employed to reduce manufacturing errors on tooth flank surfaces. Wang and Lu [24] designed a novel compounded strut with periodic structures connected in the axial direction to reduce the helicopter cabin noise generated from a gearbox.
In the existing research, the structure and mass of the shaft have always been neglected, and the modal characteristics of the system have usually been computed by considering only the gear torsional vibration; only sparse research has focused on the influence of the shaft flexibility on the vibration dynamic responses and noise of the gearbox. The stiffness and damping characteristics of shafts affect the transmissibility of vibration generated by the gears to the housing, which in turn influences the radiated noise. Therefore, it is of great significance to investigate the vibration and noise characteristics of the gearbox while considering the flexibility of the shaft under multi-source, time-varying excitation.
In this study, an improved computational method is developed to predict the vibration and noise of the gearbox. The remainder of this paper is organized as follows. Section 2 introduces the studied gearbox and the process of vibration and noise analysis. Section 3 presents the establishment of a dynamic model of a two-stage spur gear system that considers shaft flexibility. Section 4 presents the static and dynamic characteristics of the system, and the effects of the shaft flexibility on the bearing response and radiated noise of the gearbox are discussed. Sections 5 and 6 investigate the respective effects of the input speed and shaft stiffness on the vibration response and radiated noise of the gearbox. Finally, the conclusions are given in Section 7.

II. GEARBOX DESCRIPTION AND METHOD
A. THE PARAMETERS OF THE TWO-STAGE GEARBOX Fig. 1(a) and (b) are the real photos of two-stage gearbox, (a) shows the housing of gearbox, while (b) represents internal structure of the gearbox, and (c) is the schematic of the gear system. It includes two pairs of spur gears, three shafts, six deep groove ball bearings that support the shafts, and the housing. The first-stage pinion and gear are represented by p 1 and g 1 respectively. The second-stage pinion and gear are represented by p 2 and g 2 respectively. The gear parameters are listed in Table 1. Table 2 gives the shaft parameters. The ball bearings are used at the gearbox input/output shaft connections as shown in Fig 1. The input side and the output side bearings are represented by 1a, 2a, 3a, 1b, 2b and 3b respectively. Its dimensions are given in Table 3. The housing is rectangular with a height of 368 mm, width of 163 mm, and length of 365 mm. The housing wall thickness is 12 mm.   stiffness (TVMS) and gear error, the vibration is transmitted to the housing from the gears through the shafts and bearings, leading to the vibration of the housing and the acoustic emission from the housing surface. Fig. 2 presents the computational process by which to predict noise radiation from the gearbox. Multiple interconnected models are employed for the calculation of gearbox vibration transmission and acoustic radiation. The entire analysis process is composed of five components, namely the (a) gearbox dynamic excitation calculation, (b) gearbox modal analysis, (c) gearbox dynamic response analysis, (d) boundary element model construction, and (e) acoustic radiation analysis.
(a) Gearbox dynamic excitation calculation. A lumpedmass model of a spur gear pair is combined with a Timoshenko beam element model to predict the natural modes, shaft deformation, and forced response due to static transmission error excitation. The gear system is divided into many elements, including mesh elements, shaft elements, and bearing elements, for which force analysis is carried out. The entire dynamic equation is assembled by establishing different equations of motion for each element. The Newmark integration method (NIM) is utilized to obtain the bearing dynamic load. Finally, the dynamic load is exported for the calculation of the dynamic response of the gearbox.
(b) Gearbox modal analysis. A finite element model for housing is established for the modal analysis calculation, to which boundary conditions should be applied according to the assembly relation of the gearbox. The analysis outputs are the mode of vibration and the natural frequency.
(c) Gearbox dynamic response analysis. The dynamic bearing forces given by the coupled spur gear-shaft-bearing model are applied to the finite element model of the housing. To reduce the amount of computation, the vibration response and surface acceleration of the housing are computed using the mode superposition method.
(d) Boundary element model construction. The surface elements of the finite element model of the housing are extracted, and the hole on the surface of the structure is repaired by using a shell element. Finally, the acoustic boundary element model is established by remeshing the surface elements of the gearbox.
(e) Acoustic radiation analysis. To improve the calculation efficiency, an acoustic transfer vector (ATV) between the vibration of the gearbox surface and the sound field is constructed by defining the fluid attribution outside the gearbox and the sound field of the gearbox. The gearbox noise is finally computed using the boundary element model of the housing with the calculated normal velocity, and the sound pressure level at every field point is obtained.
The analysis method is subsequently described in further detail.

III. DYNAMIC MODEL OF GEARED SYSTEM
The gear transmission system is considered as an elastic system with countless degrees of freedom, where an FEM is employed. According to the shaft structure and bearing mounting positions, the continuous system is divided into a series of shaft elements, mesh elements and bearing elements, as shown in Fig 3. The shaft 1, shaft 2 and shaft 3 are respectively divided into 12 shaft elements, 8 shaft elements and 9 shaft elements. In addition, three shafts are coupled using two gear mesh elements. There are two deep groove ball bearings on each shaft to form shaft-bearing coupling elements. According to the order of shaft, the node number is defined in sequence. The first-stage pinion and gear are No.7 and No.16, the second-stage pinion and gear are No.20 and No.29. The actual determination of nodes and elements is up to the specific structure and the requirements for calculation.

A. SHAFT ELEMENT
To consider shear deformation effectively, the Timoshenko beam theory is used to analyze the shaft element. The shaft element is modeled using 2-node Timoshenko beam element with six degrees of freedom (two translation along X and Y axis and rotation along Z axis). Since this model is a parallel shaft spur gear transmission system, the excitation in the axial direction can be neglected. As shown in Fig. 4, the model of shaft element is characterized by two nodes, including i and i + 1.
The generalized displacement vector of the ith shaft element is defined as where x i , y i , x i+1 , y i+1 are translational displacements of the ith shaft element along coordinate axis, and θ zi , θ zi+1 are rotational angles of the ith shaft element around coordinate axis. The movement matrix equation of shaft element can be written as where M e is the mass matrix of the ith shaft element, K e is the stiffness matrix of ith the shaft element, C e is the damping matrix of the ith shaft element, which is defined as Rayleign damping, it can be given as where α and β are the mass and stiffness Rayleign proportional damping coefficients respectively, that can be determined by natural frequencies and modal damping ratios of the system [25].
The mass matrix and stiffness matrix are shown below The three-dimensional dynamic model of spur gear pair is shown in Fig. 5(a), where the p and g denote pinion and gear respectively. The gear pair have six degrees of freedom, which defines the coupling between the two shafts holding the gears. The pinion and gear are connected to each other by a nonlinear gear mesh spring on the plane of action (a plane tangent to both base cylinders). k m , e m , c m and α denote the time-varying mesh stiffness, mesh error, mesh damping and pressure angle of spur gear model, respectively. x p , y p , x g , y g are the translational displacements along coordinate axis of gear mesh element. θ zp and θ zg are the rotational angles around coordinate axis of gear mesh element. The generalized displacement vector of the gear mesh element is defined as The relative displacement of pinion and gear along normal line of action can be given as where e m is the mesh error of gear mesh element and given by e m = e 0 sin (2πf m1 t), V is the projective vector from generalized coordinate to normal line of action, as given in Eq.(8) where r p , r g are the base circle radius of the pinion and gear, respectively. Introducing the time-varying mesh stiffness and mesh error, the motion equation of gear mesh element can be written as where m p , m g are the mass of pinion and gear respectively, I zp , I zg are the moment of inertia of pinion and gear around the corresponding coordinate axis respectively. k m is the time-varying mesh stiffness. The solution method refers to the reference [26]. Fig. 5(b) is the solution of the time-varying stiffness and the stress cloud diagram of the single tooth meshing and double tooth meshing. In the process of gear meshing, the stress at the meshing position is the greatest, and then the stress is transmitted to the tooth root position, so the stress at the tooth root position is concentrated. The time-varying mesh stiffness curve is shown in Fig. 5(c), c m is the mesh damping of gear mesh element.
The matrix form of the motion equation of gear mesh element can be written as where M m is the mass matrix of gear element that can be written as M m = diag m p , m p , I zp , m g , m g , I zg , K m is the stiffness matrix of gear mesh element that can be written as K m = k m V T V , C m is the damping matrix of gear mesh element that can be written as C m = c m V T V , e is the mesh error vector of gear mesh element.

C. BEARING ELEMENT
The bearing supports the shaft system and transmits vibration from the gears to the housing. The bearing element is modeled using springs and dampers as shown in Fig. 6(a). The time-varying bearing stiffness is expressed in Eq. (11), and the time-varying stiffness curve as shown in Fig. 6(b).
where k a = 8.5×10 8 N /m is the static stiffness of the bearing, k 0 is the fluctuation amplitude of bearing stiffness, β 0 is the bearing phase angle, f b is the bearing passing frequency.

D. TWO-STAGE PHASE RELATIONSHIP
The phase relation of the two-stage gears should be determined before establishing a complete dynamic model of the two-stage spur gear system. A coordinate system is established as shown in Fig. 7. x p1 , y p1 , θ zp1 , x g1 , y g1 , θ zg1 are the translational displacement and torsional angle of the fist-stage pinion and gear, x p2 , y p2 , θ zp2 , x g2 , y g2 , θ zg2 are the translational displacement and torsional angle of the second-stage pinion and gear, α is the fist-stage pressure angle, β is the shaft angle among three shafts. According to the geometrical relation, the angle between the meshing line of the second-stage gear and the x p1 is reformulated as Therefore, in the second-stage gear system, the projective vector from generalized coordinate to normal line of action is given as (13) where r p2 and r g2 are the base circle radius of the second-stage pinion and gear respectively.

E. SYSTEM MODEL
When the dynamic model of each element is well defined, the overall model of the two-stage spur gear system can be determined through assembling the element dynamic models with FEM. According to the corresponding relation between the node number of each element and the node number of the gear system, the submatrix of the element is superimposed to the corresponding position of the whole matrix. The generalized displacement vector of the system is {X } = {x 1 , y 1 , θ z1 , x 2 , y 2 , θ z2 · · · x n , y n , θ zn } (n is the number of nodes, n = 1, 2, 3. . . 32). The two-stage gear transmission system is divided into 29 elements with a total of 96 degrees of freedom. The entire stiffness matrix of the two-stage parallel shaft gear system is shown in Fig. 8. The square composed of two colors is the algebraic superposition of two different types of element submatrices. The matrix form of motion equations of the system can be written as where {X (t)} is the generalized coordinate vector of the system, M , C, K are the entire mass matrix, damping matrix and stiffness matrix of the system respectively, F is the external applied force vector.

F. EXPERIMENT
The SQI wind turbine transmission system test bench is used to analyze the dynamic characteristics. As presented in Fig. 9, the time-domain signals of radial acceleration at the bearing are obtained via the acceleration sensor (356A16) installed on the end-cap of the bearing, and the data collection system produced by Dew soft is employed. The gear transmission system has parameters consistent with those of the constructed model. The input speed of the gear is varied with an increment of nearly 240 rpm, and the acceleration values are recorded under steady-state conditions. Fig. 10 presents the comparison of the measured and simulated acceleration values, from which it is evident that the measured and simulated curves agree well and the frequencies at which the resonances occur are predicted accurately. Peak A corresponds to the input speeds of 2200 rpm, where f m1 (1320 Hz) is near the third-order (1332 Hz) natural frequencies; thus, the system has strong vibration. It should be noted that the amplitude of the acceleration is found to increase significantly after passing the resonance speed area, and the acceleration changes stably with the change of the speed. However, there remain differences between the measured and simulated results because the model applies Rayleigh damping, and the vibration mode of the system changes with the change of the speed. Additionally, the coupling of air and lubricating oil also affects the system acceleration, resulting in the error between the simulation and experimental results.

IV. NUMERICAL RESULTS AND DISCUSSIONS A. NATURAL FREQUENCIES AND VIBRATION MODES
Modal analysis is a tool for the investigation of the dynamic characteristics of the system. It can be used to identify resonances. A free vibration equation of the two-stage gear system can be obtained in the following form The eigenvalues of the system can be calculated as follow where ω i , ϕ i are the natural frequencies and mode shapes of the system, respectively. Considering the effect of shaft flexibility, the natural frequencies and mode shapes of the system are calculated. Among them, the first 10 orders natural frequencies and corresponding dominant mode shapes are listed in Table 4. Fig. 11(a-c) shows the schematic of the torsional vibration of the gear, and Fig. 11(d-e) shows the schematic of the bending vibration of the shaft.
When the shaft flexibility is considered, the system has the bending vibration of the shaft. The dominant vibration of the fourth-order and eighth-order are the low-order bending mode shapes of the shaft. When the shaft flexibility is ignored, the torsional mode shape and translational mode shape of the gear pair only are analyzed. The complex mode changes of the shaft cannot be shown. Therefore, the FEM can be used for further studies.

B. SHAFT DEFORMATION ANALYSIS
Under loading conditions, the shafts are subjected to the bending moment and torque generated from the gear mesh, so the coupled bending-torsional deflection of shafts is inevitable. If the inertia force and damping force of the shaft are neglected, the force of gear shaft can be analyzed based on statics. Therefore, the deformation δ 0 can be calculated using the following equation: where K is the overall stiffness matrix of the system, F is the external applied force vector. When the torque is 100 N·m, the deformation of the shafts is shown in Fig. 12 (compared with the actual results, it is amplified by 10 5 times). The solid line represents the state VOLUME 8, 2020   after deformation and the dotted line represents the state before deformation. Results show that the maximum deformation occurs at the gear-shaft coupled node. The shaft 1 has the maximum force at the NO. 7 coupled node with the maximum deformation of 0.94 µm, shaft 2 has the maximum deformation of 1.58 µm at the NO. 20, while shaft 3 has small rotated speed and large load, so the shaft 3 has the maximum deformation of 1.92 µm at the NO. 29. Shaft deformation causes additional displacement of the gears and bearings, which deviate from the ideal position, thereby changing the meshing characteristics and system response. Therefore, it is of importance to consider shaft flexibility in gear system dynamics model.

C. BEARING DYNAMIC LOAD
Based on the two-stage gear system model proposed, dynamic responses of the system are calculated using NIM. Fig. 13 presents the schematic diagram of the bearing force analysis, and the mean value and fluctuation amplitude 85980 VOLUME 8, 2020 of each bearing dynamic load are listed in Table 5. It can be observed that the vibration of the bearings along the y-direction is significantly greater than that along the x-direction. Considering the flexibility of the shaft, due to the asymmetric layout of the shaft-gear structure, the bearing loads on both sides are unbalanced, and the dynamic load of the proximal bearing is greater than that of the distal bearing (the side near the gear is defined as the proximal side). The unbalanced bearing load results in the deflection vibration of the shafts, which can not only reduce the carrying capacity of the gears, but can also enhance the vibration and noise of the gearbox.
To illustrate the effect of the shaft flexibility on the bearing dynamic load, the bearing dynamic loads with and without regard to the shaft flexibility (SF) are calculated under the input speed of 500 rpm and torque of 100 N·m. The time-domain and frequency spectrum are shown in Fig. 14. The major composition of frequency spectrum is meshing frequency and its frequency multiplication. f m1 (300.2 Hz) and f m2 (96.7 Hz) denote the meshing frequencies of the first-stage and second-stage, respectively. f b1 (25.6 Hz) and f b2 (10.2 Hz) denote the bearing passing frequencies of bearing 1 and bearing 2, respectively. When the time-varying bearing stiffness is considered, the gear meshing frequency components and corresponding amplitude are not change at the frequency spectrum. However, there is a bearing passing frequency component at the frequency spectrum.
The bearing dynamic loads with and without regard to the shaft flexibility (SF) are compared in Fig. 14, from which it can be observed that the mean value and fluctuation amplitude of the bearing dynamic load both change. The fluctuation amplitude of the bearing dynamic load is found to be reduced after considering shaft flexibility, and the excitation frequency components and corresponding amplitude decrease in the spectrum. When the influence of the shaft flexibility is neglected, there are more frequency components, and the main frequency components of the bearing load are first 10 times that of the second-stage meshing frequency, of which 8 × f m2 has the greatest energy. When the shaft flexibility is considered, the frequency components of the bearing load on different shafts are found to be different. It is found that the main frequency components of the bearing load (bearing 1a) are f m2 , 2 × f m2 , 4 × f m1 , and 5 × f m1 , of which 2 × f m2 has the greatest energy. However, the main frequency components of the bearing load (bearing 2a) are found to be f m2 , 2×f m2 , 11×f m2 , 4×f m1 , and 5×f m1 , of which 4×f m1 has the greatest energy. Moreover, the periodic fluctuation caused by the time-varying bearing stiffness (TVBS) is weakened. The shaft flexibility therefore has substantial impacts on the bearing dynamic load that is transmitted to the housing, as it can absorb energy and decrease the vibration transmission of the gear system. The findings in this paper are consistent with those of the existing literature [6].

D. HOUSING MODAL ANALYSIS
The finite element model of the housing consist of tetrahedral elements, as shown in Fig. 15(a). The dynamic excitation of the gear system is transmitted through the bearing to the housing. In order to apply the bearing dynamic force to the housing, the mass nodes are established in the center of the bearing holes. Then, the coupling relation between the mass nodes and the inner surface nodes of the bearing holes are defined to support the inner holes of bearing. At the same time, the bottom of the housing is fully restrained. Finally, modal analysis of housing is carried out.
The first five order natural frequencies and vibration modes of housing are shown in Fig. 15(b-f). Because the lower part of housing is constrained, its vibration amplitude is smaller, while the vibration on the upper part of housing is more intense. The first-order vibration mode is bending on the top and swaying on front and rear end of housing; the second-order vibration mode is swaying on left and right sides of housing; the third-order vibration mode shows torsional vibration mode; the fourth-order and the fifth-order vibration mode are axial vibration on bearing pedestal.

E. BE MODEL AND ATV
According to the correlation between the acoustic boundary element model and the finite element model of the housing,   the bearing hole on the surface of the housing is repair by using the shell element, and the acoustic boundary element model of the housing is established by remeshing the surface elements of the housing.
The hemispherical sound field is established with the gearbox as the center and 1 m as the radius, as shown in Fig. 16(a).
To improve efficiency and reduce unnecessary calculation, acoustic transfer vector (ATV) is utilized to establish the corresponding relationship between the vibration velocity on the surface of the housing and the sound pressure at the field point, as shown in Fig. 16(b). According to the boundary element model of housing and sound field, the ATV is calculated. For calculation, the sound speed is 341 m/s, the air density is 1.21 kg/m3, and the step pitch is 10 Hz. The ATV is calculated only once and then used in conjunction with separately calculated bearing load.
{v n } is the surface velocity vector, P is the sound pressure level of field point where ATV is the acoustic transfer vector, ω is angular frequency.

F. RADIATED NOISE OF GEARBOX
The fluctuation amplitude of the bearing dynamic load (removing the mean value) is applied to the mass node in the center of the bearing hole to calculate the vibration response of the gearbox. The velocity of housing surface is considered to be the boundary condition, the radiated noise of the gearbox is calculated with ATV. The noise spectrum is calculated under the input speed of 500 rpm and torque of 100 N·m, as presented in Fig. 17(a). It is observed from the noise distribution that the noise on the top field point much lower than that on the right and left sides. Meanwhile, the field point on both sides of the bearing are in symmetric distribution, so the sound pressure level of both sides are equal. From the sound pressure level curves, it is well noticed that peaks are generated on the positions of f n1 , 3×f m2 , 6×f m2 , 9×f m2 , 4×f m1 , 5×f m1 , and 6×f m1 (f n1 = 129.7 Hz denotes the first-order natural frequency of the gear system, f m1 = 300.2 Hz and f m2 = 96.7 Hz denote the meshing frequencies of the first-stage and second-stage gear pair, respectively). The corresponding frequency of the maximum radiated noise is 4 × f m1 (1201 Hz), which is close to the second-order natural frequency of housing (1204.1 Hz). From modal analysis of housing, the second-order mode shape is concaving and swaying on left and right sides of housing. Therefore, the noise on both sides of the gearbox is significantly larger than that on the top.
The shaft flexibility affects the transmissibility of the vibration and has a significant effect on the bearing dynamic load. To illustrate the effect of the shaft flexibility on the radiated noise of the gearbox, the noise (both with and without SF) on the top field point and right field point are compared, as presented in Figs. 17(b-c). The results demonstrate that, when the flexibility of the shaft is considered, the radiated noise of the gearbox increases. Although the flexibility of the shaft can absorb energy and decrease the transmission of the vibration of the gear system, the unbalanced load of the bearing at both ends of the shaft causes the shaft to undergo deflection vibration, which increases the vibration and ultimately results in the increase of the radiated noise of the gearbox. When the flexibility of the shaft is ignored, the sound pressure levels around 200-900 Hz and near 1800 Hz notably increase, which can be attributed to the bearing having more excitation frequency components and a larger vibration energy. In addition, when the excitation frequency component is close to the structural natural frequency, the vibration amplitude and the sound pressure level are found to increase. Therefore, 9 × f m2 (870 Hz) is closer to the first-order natural frequency of the housing (863.7 Hz), resulting in the greater vibration of the gearbox and greater noise radiation. In the first-order vibration mode, the top of the housing vibrates violently, so the noise of the top field point is greater than that of both sides.

V. EFFECT OF INPUT SPEED
The change of input speed not only influences gear meshing condition, but also give rise to the change of gear meshing frequency and harmonics. In this section, in order to investigate the effects of input speed on dynamic respond and radiated noise of gearbox, the bearing load and radiated noise at input speeds from 500 rpm to 6000 rpm are calculated.

A. BEARING LOAD
The bearing dynamic loads at different input speeds are calculated using the proposed model. The comparison of bearing load fluctuation amplitude (BLFA) at different input speeds is presented in Fig. 18. It can be determined that the BLFA has a significant increasing trend with the increase of input speed, and many peaks are clearly visible. In addition, bearing load at both ends of gear has the same changing trend. The bearing 1a and bearing 1b vibrate strongly at the input speeds of 2000 rpm, 3100 rpm, 4100 rpm, 4700 rpm, and 4900 rpm, and the main excitation frequency components at these peaks are listed in Table 6. With the increase of input speed, the amplitude of the high-frequency components gradually decrease, and the low-frequency components gradually become the main excitation frequency component of the system. The result is completely consistent with the finding in [27]. When the input speed is 2000 rpm, 3100 rpm and 4100 rpm, the main excitation frequency components (6×f m2 , 4×f m2 , 3×f m2 ) are close to the seventh-order natural frequency of the gear system (2390 Hz), and the system has greater vibration energy. The vibration mode is the torsional and translational vibration of the second gear pair. When the input speed is 4700 rpm and 4900 rpm, due to the bending vibration mode of the shaft 1 (2883 Hz) is sensitive to the excitation, the vibration of the bearing at both ends of the gear is intense, and the bearing load difference also increase suddenly. The results demonstrate that the bending vibration mode of the shaft is added to the system, and it is conducive to accurately reflect the dynamic characteristics of the gear system.

B. RADIATED NOISE OF GEARBOX
To quantify the effect of the input speed on the radiated noise of the gearbox, the effective sound pressure level (ESPL) is calculated to show the noise levels at different input speeds, as shown in Eq. (19).
where P(t) is the sound pressure, P 0 is the reference sound pressure, T is the solution period.
The main excitation frequency components of gearbox are composed of gear mesh frequency and its harmonics, and the vibration amplitude is affected by the natural frequency of gear transmission system and housing.
Noise spectrum of the gearbox in the input speed range from 500 rpm to 6000 rpm is presented in Fig. 19. It can be clearly determined that the noise curves of each field point exhibit a gradually increasing trend, and the noise at the top field point is gradually larger than that at the two sides. This is because the vibration amplitude of the bearing gradually increase with the increase of the input speed, and the deflection vibration of the shaft becomes more and more violent; thus, the vibration of the gearbox increase significantly at the top field point, leading to an increase in the radiated noise of the gearbox. When the input speeds is 500 rpm and 900 rpm, the main excitation frequency components are 4 × f m1 and 7 × f m2 , respectively, which are closer to the second-order natural frequency of the housing (1204 Hz). The corresponding vibration mode is concave vibration on the left and right sides. Therefore, the noise on the top field point is much lower than that on the right and left sides. At the input speeds of 2000 rpm, 3100 rpm and 4100 rpm, the bearing dynamic load is found to fluctuate violently, and the main frequency components are 6 × f m2 , 4 × f m2 , and 3 × f m2 , respectively, which are closer to the fifth-order natural frequency of the housing (2298.7 Hz). Moreover, the vibration amplitude of the gearbox increases due to the deflection vibration of the shaft; therefore, the noise on the top of the gearbox is significantly larger than that on both sides. Although the bending vibration mode of the input shaft is found to cause dramatic fluctuations in the bearing load at the input speeds of 4700 rpm and 4900 rpm, the natural mode of the gearbox is not sensitive to the excitation, so the radiated noise of the gearbox does not increase. At the input speed of 5400 rpm, because the eighth-order natural frequency of gearbox (3101 Hz) is highly sensitive to the excitation, the vibration of the gearbox mainly occurs at the top; thus, the noise on the top field point is much larger than that on both sides.

VI. EFFECT OF SHAFT STIFFNESS
The shaft flexibility affects the transmissibility of the vibration and the radiated noise of the gearbox. In the design of the low-noise gearbox, it is vitally essential to explore the effects of the shaft stiffness on dynamic response and radiated noise of gearbox. The reasonable design of the shaft stiffness can avoid the system resonance and reduce the vibration and noise of the gearbox. The shaft stiffness be changed by modifying its diameter or shear modulus. In this section, the effects of the shaft stiffness on static and dynamic behaviors of the two-stage gearbox are investigated.

A. NATURAL FREQUENCIES AND VIBRATION MODES
The natural frequencies and dominate vibration modes of the gear transmission system with shaft stiffness of K , 2 × K , 3 × K , 4 × K , 5 × K , 6 × K , and 7 × K (initial stiffness is K ) are calculated respectively. The natural frequencies of the gear transmission system are greatly affected by the change of shaft stiffness, while corresponding vibration modes are less affected. The first 6 orders natural frequencies of the system are presented in the Fig. 20, from which it can be seen that the natural frequencies of the system increases monotonically with the increase of the shaft stiffness. This is because the stiffness of the system increases with the increase of the shaft stiffness. When the shaft stiffness reaches 3×K , the third-order dominate vibration modes of the system change from torsional vibration of the gear to the bending vibration of the shaft. For the higher-order vibration modes of the system (16th, 23th, 24th, 30th, 31th), the bending vibration of the shaft is transformed into torsional vibration of the gear, as a result, the number of the torsional vibration modes and vibration energy of the gear transmission system increase.

B. SHAFT DEFORMATION
The static deformation of the shaft 1, shaft 2, and shaft 3 with the shaft stiffness of K , 2 × K , 3 × K , 4 × K , 5 × K , 6 × K and 7 × K are calculated respectively. As presented in the Fig. 21, the deformation of gear-shaft coupled nodes are defined as X mp1 , X mg1 , X mp2 , and X mg2 , respectively. It can be found that the deformation of the shaft decreases suddenly followed by a slow decrease, each node shows the same trend. When the shaft stiffness increases from K to 2×K , X mp1 , X mg1 , X mp2 , and X mg2 decrease by 44%, 41.7%, 43.7% and 44.1%, respectively. The results show that gear meshing characteristics are less and less affected by shaft deformation with the increase of the shaft stiffness.

C. BEARING LOAD
The change of the shaft stiffness affects the deformation of the shaft, and the reduction of the shaft deformation affects the distribution of the bearings load at both ends. Therefore, the bearing loads of the two-stage gear system with the different shaft stiffness are calculated using the proposed model. The bearing load difference between the two ends of the gear is defined as follows F y = F ay − F by (20) where F ay and F by denote the mean load of bearing a and bearing b in the y direction, respectively. The bearing load difference with different shaft stiffness is given in Fig. 22, from which it can be observed that with the increase of the shaft stiffness, the loads of bearing a and bearing b gradually become equal. This can be attributed to the gradual reduction of the deformation of the shaft give rise to the dynamic loads of the bearing 1a, 2b and 3b gradually decrease, while the dynamic loads of the bearing 1b, 2a and 3a gradually increase. It means that the deflection vibration of the shaft and noise of gearbox gradually reduce with the increase of the shaft stiffness. Fig. 23(a) presents the curves of the bearing load fluctuation amplitude (BLFA) at different shaft stiffness. With the multiple increases of the shaft stiffness, the BLFA curves exhibit a slow increasing trend with significant fluctuation. This can be attributed to the gradual weakening of the capacity of the shaft to absorb the vibration energy of the system with the increase of the stiffness of the shaft, and the increasing excitation frequency components and amplitudes of the bearing load. It can be clearly found that when the stiffness of the shaft is K , the vibration of the bearings is small. When the shaft stiffness reaches 3 × K and 6 × K , the bearings vibrate violently, leading to the stronger vibration of the system. To clarify the origins of this phenomenon, the waterfall of the load of bearing2a at different shaft stiffness is presented in Fig. 23(b). With the increase of the shaft VOLUME 8, 2020  stiffness, increasingly more multiple-frequency components occur, and the amplitude increases. From the preceding analysis, it can be determined that that the increase of the shaft stiffness leads to the increase of the natural frequency of the system and the change of the vibration mode. When the shaft stiffness reaches 6×K , the main frequency components are first 7 times that of the second-stage meshing frequency in the frequency spectrum. In addition, the change of the shaft stiffness results in the changes in the amplitudes of the excitation frequency components. When the shaft stiffness reaches 3 × K and 6 × K , the main excitation frequency components are 5 × f m1 (1501 Hz) and 2 × f m1 (600 Hz), respectively. This is because they are respectively close to the third-order (1513 Hz) and first-order (608 Hz) natural frequencies of the gear system, leading to the resonance of the system. Therefore, the increase of the shaft stiffness gives rise to stronger gearbox vibration, which affects the radiated noise of the gearbox.

D. RADIATED NOISE
The effective sound pressure levels (ESPLs) of the field points on the top and both sides of the housing are calculated at different shaft stiffness, as presented in Fig. 24. The noise of the field points on both sides and the top of the gearbox exhibits a nonlinear change. In addition, the radiated noise of the housing changed from one pattern (the maximum noise is at both sides of the housing) to another pattern (the maximum noise is at the top of the housing). Although the fluctuation amplitude of the bearing load increases with the increase of the shaft stiffness, the radiated noise of the gearbox decreases with the increase of the shaft stiffness. This is because the deflection vibration of the shafts is weakened and the dominate resonance modes of the housing is changed.
The noise spectra at different shaft stiffness are presented in Fig. 25, from which it is evident that the noise distribution changes due to the change of the shaft stiffness. The distribution curves of the sound pressure level are found to peak at 3 × f m2 , 2 × f m1 , 3 × f m1 , 4 × f m1 , 5 × f m1 , and 6 × f m1 . When the shaft stiffness is K , more intense noise  is found to have appear at 4 × f m1 (1201 Hz). Because it has the maximum vibration energy and is consistent with the second-order natural frequency of housing (1204 Hz), the sound pressure level at this frequency is larger than that of other stiffness. When the shaft stiffness reaches 3 × K , the excitation frequency component (1501 Hz) is found to have the highest noise because it is closer to the third-order natural frequency of the gear transmission system (1513 Hz), and its vibration shape is the bending vibration of shaft 3. Therefore, the noise at the right field point is greater than that at the left field point. When the shaft stiffness is 6 × K , intense noise is found to appear at around 600 Hz and 900 Hz due to the excitation frequencies are close to the first-order natural frequency of the gear system (608 Hz) and housing (864 Hz), resulting in an increase in the vibration and noise of the gearbox.

VII. CONCLUSION
In this paper, a coupled dynamic model of a spur gearshaft-bearing system that considers the flexibility of the shaft is developed based on the Timoshenko beam theory. A novel computational method for a two-stage gearbox is proposed to address the vibro-acoustic propagation of the gearbox. The influences of the shaft flexibility, input speed, and shaft stiffness on the vibration response and radiated noise of the two-stage gearbox are investigated. The main conclusions of this research can be summarized as follows.
(1) The dynamic model of the gear system is validated with experiments. By modeling the mass and structure of the shafts, the vibration modes of the gear transmission system are found to undergo the bending vibration of the shaft. The deformation of the shaft is found to result in the difference of the bearing load at both ends of the gear, which leads to the deflection vibration of the shaft.
(2) When considering the flexibility of the shaft, the fluctuation amplitude and excitation frequency components of the bearing dynamic load are found to decrease significantly, and the radiated noise is enhanced. The main resonance mode of the gearbox is found to change, and the noise near the low-frequency band (200-900 Hz) and the high-frequency band (1800 Hz) is notably reduce.
(3) As the input speed increases, the low-frequency component gradually becomes the main excitation frequency component of the system, and the increasing dynamic load results in the stronger vibration and noise of the gearbox. With the increase of the input speed, the bending vibration mode of the shaft is added to the system, resulting in the increase of the vibration amplitude. Additionally, the noise at the top of the gearbox gradually becomes much larger than that of the two sides.
(4) Shaft stiffness has a significant effect on the vibration response and radiated noise of the gearbox. With the increase of the shaft stiffness, the natural frequency of the system increases and the shaft deformation decreases, and the bearing loads at both ends of the gear gradually become equal. The fluctuation amplitude of the bearing dynamic load increases, but the radiated noise of the gearbox decreases due to the weakening of the deflection vibration of the shaft.