Quasi-Static Characteristics and Vibration Responses Analysis of Helical Geared Rotor System with Random Cumulative Pitch Deviations

: As one of the long period gear errors, the e ﬀ ects of random cumulative pitch deviations on mesh excitations and vibration responses of a helical geared rotor system (HGRS) are investigated. The long-period mesh sti ﬀ ness (LPMS), static transmission error (STE), as well as composite mesh error (CMS), and load distributions of helical gears are calculated using an enhanced loaded tooth contact analysis (LTCA) model. A dynamic model with multi degrees of freedom (DOF) is employed to predict the vibration responses of HGRS. Mesh excitations and vibration responses analysis of unmodiﬁed HGRS are conducted in consideration of random cumulative pitch deviations. The results indicate that random cumulative pitch deviations have signiﬁcant e ﬀ ects on mesh excitations and vibration responses of HGRS. The curve shapes of STE and CMS become irregular when the random characteristic of cumulative pitch deviations is considered, and the appearance of partial contact loss in some mesh cycles leads to decreased LPMS when load torque is relatively low. Vibration modulation phenomenon can be observed in dynamic responses of HGRS. In relatively light load conditions, the amplitudes of sideband frequencies become larger than that of mesh frequency and its harmonics (MFIHs) because of relatively high contact ratio. The inﬂuences of random cumulative pitch deviations on the vibration responses of modiﬁed HGRS are also discussed.


Introduction
Due to high power density, high reliability, as well as high transmission efficiency, gear transmission devices are widely applied to wind turbine, marine, automotive, and other industrial fields. Owing to manufacturing and heat treatment processes, which cause tooth displacement relative to its ideal position, manufacturing errors with various distribution forms and different magnitudes can be found in actual gears.
During the past several decades, a great deal of researchers have investigated the dynamic behaviors of spur/helical gear systems (in consideration of manufacturing errors) in published works. Gear errors consist of short and long period components. The short period gear errors were usually assumed as simple harmonic functions related to mesh frequency, the long period ones were regarded as a sine wave related to shaft frequency, and the time-varying mesh stiffness was regarded as a stiffness excitation which was independent with gear errors [1][2][3][4][5]. Moreover, when gear errors are introduced, the three-dimensional contact state of the tooth surface will differ under different load conditions. The phenomenon of partial contact loss may occur and lead to the changes of mesh stiffness and transmission error [6][7][8][9]. Hence, there is a nonlinear coupling relationship between mesh stiffness and gear errors because the contact patterns are dependent with load conditions. The short period gear errors mainly include profile deviation, helix deviation, and gear modification. Mesh misalignment [10], which is caused by mounting errors and system deflection, is also regarded as a special short period gear error. Some scholars established the calculation models of load distribution [11][12][13][14], mesh stiffness, and transmission error [15][16][17][18] in consideration of short period gear errors, and further numerically and experimentally investigated the dynamic behaviors of gear system [19][20][21].
Cumulative pitch deviation is a very important and inevitable long period gear error. When introducing the long period gear error, the gear engagement processes will differ during different mesh cycles, and the system will be subjected to hybrid mesh excitations related to the mesh frequency and shaft frequency. In published works, considering cumulative pitch deviations, the determination of long period mesh excitations using three-dimensional contact analysis and dynamic behaviors analysis of gear system are rather rare. In recent years, few relevant research was conducted, and most of the research focused on the spur gear pair. Fernández [22] developed an enhanced model to study the influence of cumulative pitch deviation on load sharing ratio, as well as static transmission error (STE) and dynamic mesh forces of a spur gear pair. The results showed that index errors not only influence the vibration and noise, but also lead to high overloads. Talbot [23] and Inalpolat [24] employed a special test setup with transmission ratio of 1:1 to study the effect of tightly controlled intentional cumulative pitch deviations on the baseline vibration responses and vibration frequency spectrum of spur gear pairs. Moreover, Yuan [25] preliminarily investigated the influences of sinusoidal cumulative pitch deviations on mesh excitations and frequency spectrum of a helical gear pair.
In engineering practices, the pitch deviation values distributed over gear teeth usually has random characteristic because of manufacturing tolerance. For helical gear pairs, the number of teeth engaged at the same time is more than spur gear pairs. In addition, gear tooth flank is usually modified in order to reduce vibration. At this time, gear modification and cumulative pitch deviation have a comprehensive influence on tooth engagement process of helical gear pairs. Thus, the effects of cumulative pitch deviations on long period mesh excitations and vibration responses of the helical geared rotor system (HGRS) are much more complex. The relevant researches have not been reported in published works.
As the above analysis of the published research, the influence mechanism of cumulative pitch deviations on long period mesh excitations and vibration responses of the unmodified and modified HGRS is still unclear. It is necessary and important to investigate the changing law of long-period mesh stiffness (LPMS), STE, and composite mesh error (CMS) under the influence of random cumulative pitch deviations, and deeply explore the vibration spectrum characteristics of the unmodified and modified HGRS through introducing these long period mesh excitations, which are determined using three-dimensional contact analysis. The rest of this paper is organized as follows: Section 2 gives an enhanced loaded tooth contact analysis (LTCA) model of helical gear pairs with cumulative pitch deviations for determining LPMS, STE, and CMS, and establishes a multi degrees of freedom (DOF) dynamic model for HGRS using Timoshenko beam theory and finite element method. Section 3 gives some numerical results of mesh excitations and vibration responses of unmodified HGRS with random cumulative pitch deviations, and further analyzes the influence of random cumulative pitch deviations on the dynamic responses of modified HGRS. Section 4 draws some conclusions from the numerical results. The detailed flow chart of the study is given in Figure 1.

Description of Cumulative Pitch Deviations for Helical Gears
It is well known that cumulative pitch deviation and pitch deviation co-exist in actual gears. Pitch deviation denotes the circumferential position error of one tooth flank relative to the previous one. Cumulative pitch deviation refers to the algebraic difference between actual and theoretical arc length of any teeth spaces. Due to cumulative pitch deviations, the actual load sharing ratio between gear teeth is quite different from theoretical one. The phenomenon of overload or partial contact loss may occur in some potential contact region, and the actual contact ratio of the gear pair will be lower or higher than theoretical one.
Cumulative pitch deviation and pitch deviation are both measured along the tangential direction of standard pitch circle of gears. Figure 2a shows a helical gear with cumulative pitch deviations, and Figure 2b gives the mating state of the gear pair with cumulative pitch deviation when the contact tooth number is five in one engagement position. Here, the gear teeth numbered as zero is the reference tooth, bt P and pt f are theoretical tooth pitch and pitch deviation, respectively. ′ btk P and pk F are actual tooth pitch and cumulative pitch deviation for the gear tooth numbered as k, respectively. The pitch deviation ptk f of the k-th gear tooth is equal to the difference of two adjacent cumulative pitch deviations, it can be written as In three-dimensional contact analysis model, the initial normal clearances between potential contact point pairs denote gear errors. The error values of cumulative pitch deviations and pitch deviations must be converted to the normal direction of tooth surface as follows cos cos where pbn f is the pitch deviation along the normal direction of gear tooth flank, α t is the transverse pressure angle, β b is the helix angle of helical gears.

Description of Cumulative Pitch Deviations for Helical Gears
It is well known that cumulative pitch deviation and pitch deviation co-exist in actual gears. Pitch deviation denotes the circumferential position error of one tooth flank relative to the previous one. Cumulative pitch deviation refers to the algebraic difference between actual and theoretical arc length of any teeth spaces. Due to cumulative pitch deviations, the actual load sharing ratio between gear teeth is quite different from theoretical one. The phenomenon of overload or partial contact loss may occur in some potential contact region, and the actual contact ratio of the gear pair will be lower or higher than theoretical one.
Cumulative pitch deviation and pitch deviation are both measured along the tangential direction of standard pitch circle of gears. Figure 2a shows a helical gear with cumulative pitch deviations, and Figure 2b gives the mating state of the gear pair with cumulative pitch deviation when the contact tooth number is five in one engagement position. Here, the gear teeth numbered as zero is the reference tooth, P bt and f pt are theoretical tooth pitch and pitch deviation, respectively. P btk and F pk are actual tooth pitch and cumulative pitch deviation for the gear tooth numbered as k, respectively. The pitch deviation f ptk of the k-th gear tooth is equal to the difference of two adjacent cumulative pitch deviations, it can be written as

Determination of LPMS, STE, and CMS
In our previous work, a basic LTCA model with high computational efficiency has been developed [26]. It can be used to calculate the short period mesh stiffness and transmission error of a cylindrical gear pair with short period tooth flank errors. In this section, considering the difference of tooth engagement process in different mesh cycles and the actual engagement condition of each tooth pair, the basic LTCA model is extended to an enhanced LTCA model, which can be applied to the calculation of long period mesh excitations. Figure 3 shows the enhanced LTCA model for a helical gear pair during one mesh cycle, and gives the arrangement of potential contact lines and contact points on plane of action (POA) B1B2B3B4. For a helical gear pair with total contact ratio of 2 to 3, it can be divided into 2 teeth and 3 teeth zone. rp and rg are the radius of base circle of driving and driven gear. Op and Og are the rotational center of the two gears. βb is helix angle of gear base circle. The rotational speeds of the two gears are denoted as ωp and ωg. The Cartesian coordinate system XMOMYM is defined for LTCA calculation. Figure 4 shows the initial contact condition of one contact line pair before load for unmodified and modified gear pairs in consideration of cumulative pitch deviations. fpbn1 and fpbn2 are the pitch deviations for driving and driven gear respectively, Mi1 and Mi2 are the modification values for contact point i on driving and driven gear.   In three-dimensional contact analysis model, the initial normal clearances between potential contact point pairs denote gear errors. The error values of cumulative pitch deviations and pitch deviations must be converted to the normal direction of tooth surface as follows where f pbn is the pitch deviation along the normal direction of gear tooth flank,α t is the transverse pressure angle, β b is the helix angle of helical gears.

Determination of LPMS, STE, and CMS
In our previous work, a basic LTCA model with high computational efficiency has been developed [26]. It can be used to calculate the short period mesh stiffness and transmission error of a cylindrical gear pair with short period tooth flank errors. In this section, considering the difference of tooth engagement process in different mesh cycles and the actual engagement condition of each tooth pair, the basic LTCA model is extended to an enhanced LTCA model, which can be applied to the calculation of long period mesh excitations. Figure 3 shows the enhanced LTCA model for a helical gear pair during one mesh cycle, and gives the arrangement of potential contact lines and contact points on plane of action (POA) For a helical gear pair with total contact ratio of 2 to 3, it can be divided into 2 teeth and 3 teeth zone. r p and r g are the radius of base circle of driving and driven gear. O p and O g are the rotational center of the two gears. β b is helix angle of gear base circle. The rotational speeds of the two gears are denoted as ω p and ω g . The Cartesian coordinate system X M O M Y M is defined for LTCA calculation. Figure 4 shows the initial contact condition of one contact line pair before load for unmodified and modified gear pairs in consideration of cumulative pitch deviations. f pbn1 and f pbn2 are the pitch deviations for driving and driven gear respectively, M i1 and M i2 are the modification values for contact point i on driving and driven gear. the two gears. βb is helix angle of gear base circle. The rotational speeds of the two gears are denoted as ωp and ωg. The Cartesian coordinate system XMOMYM is defined for LTCA calculation. Figure 4 shows the initial contact condition of one contact line pair before load for unmodified and modified gear pairs in consideration of cumulative pitch deviations. fpbn1 and fpbn2 are the pitch deviations for driving and driven gear respectively, Mi1 and Mi2 are the modification values for contact point i on driving and driven gear.  In the LTCA model, the contact problem of gear tooth flanks in each engagement position is regarded as a line contact problem of two elastic bodies. Each contact line is evenly discretized into several contact points so that the line contact problem will be transformed into point contact problem. The gear errors on tooth flank will be introduced into potential contact point pairs as initial contact clearances. When the load torque is applied, the tooth flank of mating gear teeth will come into contact, and the gear deflections will occur. Gear deflections consists of global deflections of mating gear teeth and local contact deformation of contact points. The global deflection changes linearly with the load, it can be determined using the flexibility matrix of contact points and loads on these contact points. The finite element substructure method can be employed to determine the global deflection flexibility of contact points [26]. The contact deformation of each contact point changes nonlinearly with the load, which can be calculated using analytical formula of contact deformation [15]. Introducing cumulative pitch deviations and tooth surface modifications into the initial clearances between contact point pairs, the contact condition of mating gear teeth in each engagement position can be expressed as [15] δ χ ξ ω For each engagement position, the nonlinear matrix Equation (3) can be solved using an iteration contact algorithm, which contains two control loops [26]. When the solutions in each engagement position of all mesh cycles are solved one by one, the long-period load distribution and STE will be determined.
The LPMS of helical gear pairs with cumulative pitch deviations can be determined using [15] N N F Figure 4. Initial contact condition before load in consideration of cumulative pitch deviations.
In the LTCA model, the contact problem of gear tooth flanks in each engagement position is regarded as a line contact problem of two elastic bodies. Each contact line is evenly discretized into several contact points so that the line contact problem will be transformed into point contact problem. The gear errors on tooth flank will be introduced into potential contact point pairs as initial contact clearances. When the load torque is applied, the tooth flank of mating gear teeth will come into contact, and the gear deflections will occur.
Gear deflections consists of global deflections of mating gear teeth and local contact deformation of contact points. The global deflection changes linearly with the load, it can be determined using the flexibility matrix of contact points and loads on these contact points. The finite element substructure method can be employed to determine the global deflection flexibility of contact points [26]. The contact deformation of each contact point changes nonlinearly with the load, which can be calculated using analytical formula of contact deformation [15]. Introducing cumulative pitch deviations and tooth surface modifications into the initial clearances between contact point pairs, the contact condition of mating gear teeth in each engagement position can be expressed as [15] Appl. Sci. 2020, 10, 4403 6 of 18 where [δ G ] denotes the summation of flexibility matrix of the mating gear teeth, {F} refers to the load distribution vector on contact points after load, {χ C } is the contact deformations of contact points, {ω} is the residual clearances between contact point pairs after load, {ξ} is the initial clearances between contact point pairs before load. P is the normal force of gear pair.
For each engagement position, the nonlinear matrix Equation (3) can be solved using an iteration contact algorithm, which contains two control loops [26]. When the solutions in each engagement position of all mesh cycles are solved one by one, the long-period load distribution and STE will be determined.
The LPMS of helical gear pairs with cumulative pitch deviations can be determined using [15] The CMS of helical gear pairs can be calculated using [15]

Dynamic Model of HGRS
The research in this paper is aimed at cumulative pitch deviations, which have long periods, and the motivation is to study the influence mechanism of cumulative pitch deviations on the mesh excitations and vibration spectrum characteristics of HGRS. Thus, it is unnecessary to select complex gear transmission device, and a single stage gear system can satisfy the requirement. In the published researches, a single-stage gear transmission device with transmission ratio of 1:1 is used for several times [23,24,[27][28][29]. In order to avoid too large s computational effort, according to those published works, a HGRS with transmission ratio of 1:1 is selected as the research object. Figure 5 shows the structure diagram of a HGRS. Shaft 1 and shaft 2 are solid shafts where the driving and driven gear are located. Both shafts are composed of several shaft segments with different diameters and lengths. The parameters of the helical gear pair are given in Table 1. Each shaft is composed of 5 shaft segments. The outer diameters of shaft segments of the two shafts from left to right are 70 mm, 90 mm, 110 mm, 90 mm, and 70 mm, respectively. For shaft 1, the length of shaft segments from left to right is 100 mm, 110 mm, 220 mm, 110 mm, and 60 mm, respectively. For shaft 2, the length of shaft segments from left to right is 60 mm, 110 mm, 220 mm, 110 mm, and 100 mm, respectively. Considering shaft structure, gear and bearing positions, each shaft is divided into 13 elements, as shown in Figure 5. Three basic elements (gear mesh element, shaft element, and bearing element) are used to build the dynamic model of HGRS.

Dynamic Model of HGRS
The research in this paper is aimed at cumulative pitch deviations, which have long periods, and the motivation is to study the influence mechanism of cumulative pitch deviations on the mesh excitations and vibration spectrum characteristics of HGRS. Thus, it is unnecessary to select complex gear transmission device, and a single stage gear system can satisfy the requirement. In the published researches, a single-stage gear transmission device with transmission ratio of 1:1 is used for several times [23,24,[27][28][29]. In order to avoid too large s computational effort, according to those published works, a HGRS with transmission ratio of 1:1 is selected as the research object. Figure 5 shows the structure diagram of a HGRS. Shaft 1 and shaft 2 are solid shafts where the driving and driven gear are located. Both shafts are composed of several shaft segments with different diameters and lengths. The parameters of the helical gear pair are given in Table 1. Each shaft is composed of 5 shaft segments. The outer diameters of shaft segments of the two shafts from left to right are 70 mm, 90 mm, 110 mm, 90 mm, and 70 mm, respectively. For shaft 1, the length of shaft segments from left to right is 100 mm, 110 mm, 220 mm, 110 mm, and 60 mm, respectively. For shaft 2, the length of shaft segments from left to right is 60 mm, 110 mm, 220 mm, 110 mm, and 100 mm, respectively. Considering shaft structure, gear and bearing positions, each shaft is divided into 13 elements, as shown in Figure 5. Three basic elements (gear mesh element, shaft element, and bearing element) are used to build the dynamic model of HGRS.    Introducing spring and clearance, gear mesh element is used to simulate instantaneous engagement state of gear pairs. It consists of two nodes, each node has 6 DOF, as shown in Figure 6. The spring stiffness k m refers to the LPMS, and the clearance e m denotes the CMS. The two gears rotate around O p and O g . β b and α are the base helix angle and mesh angle. ψ is the phase angle.
mesh element can be written as  Timoshenko beam element is used to establish the dynamic model of shaft element, which includes 2 nodes with 12 DOF [15]. Ignoring the housing flexibility, bearing element is established using bearing stiffness and damping which connects shaft with foundation. In general, owing to the number change of loaded rollers, the stiffness of rolling bearing is usually time-dependent, whereas the published work shows that the numerical simulation results of a geared rotor system agree well with the experimental results when constant bearing stiffness is used. Thus, the time-dependent characteristic of the rolling bearing is neglected in this study. The stiffness matrix, damping matrix, and motion equation of bearing element can be found in [15]. The generalized coordinate vector of gear mesh element can be defined as u m = x p , y p , z p , θ xp , θ yp , θ zp , x g , y g , z g , θ xg , θ yg , θ zg T The relative displacement of driving and driven gear along normal line of action can be calculated using τ m = Tu m where T can be expressed as Appl. Sci. 2020, 10, 4403 8 of 18 Based on the Newton motion law, introducing LPMS and CMS, the motion equations for gear mesh element can be written as The matrix form of Equation (8) can be expressed as where M m and e m (t) are mass matrix and CMS vector, respectively. Timoshenko beam element is used to establish the dynamic model of shaft element, which includes 2 nodes with 12 DOF [15]. Ignoring the housing flexibility, bearing element is established using bearing stiffness and damping which connects shaft with foundation. In general, owing to the number change of loaded rollers, the stiffness of rolling bearing is usually time-dependent, whereas the published work shows that the numerical simulation results of a geared rotor system agree well with the experimental results when constant bearing stiffness is used. Thus, the time-dependent characteristic of the rolling bearing is neglected in this study. The stiffness matrix, damping matrix, and motion equation of bearing element can be found in [15].

System Model
When the dynamic model for each element is established, the finite element method is employed to assemble the element models according to the system structural parameters, and the dynamic model of the HGRS can be obtained, as shown in Figure 7.

System Model
When the dynamic model for each element is established, the finite element method is employed to assemble the element models according to the system structural parameters, and the dynamic model of the HGRS can be obtained, as shown in Figure 7. The matrix form of motion equations of the HGRS can be expressed as where G q is the generalized coordinates of HGRS, where q G is the generalized coordinates of HGRS, M G , C G , K G are mass matrix, damping matrix and stiffness matrix of HGRS, respectively, e G denotes the CMS vector of HGRS, F G refers to load vector of HGRS.
In the actual operating condition, the gear pair is subjected to dynamic mesh force. For a spur gear system, it is easy to observe the strong nonlinear dynamic behavior ("jump" phenomenon), which is caused by tooth separations under resonant condition. However, for a helical gear system, the fluctuation of dynamic mesh force is far less than that of spur gear pairs owing to the much higher contact ratio. The published experimental work shows that strong nonlinear phenomenon ("jump" phenomenon) caused by tooth separations usually does not occur even under resonant condition for a helical gear system [21,30]. Moreover, it is suggested that the dynamic responses of a helical gear system, which is calculated by using quasi-static mesh excitations, is very close to the experimental results in the published works [15,30]. The research in this paper focuses on cumulative pitch deviations, which have a long period, and the established dynamic model of HGRS has large DOF. Considering the huge computational effort of the three-dimensional dynamic contact model and the conclusions in the published works [15,21,30], the quasi-static mesh excitations obtained from the enhanced LTCA model are used to predict the vibration responses of the HGRS.

Numerical Results and Discussion
In engineering practices, the cumulative pitch deviations with sine wave, which are used in published works [22,31], rarely occurs because of manufacturing tolerances. It is more practical to have random cumulative pitch deviations on gear teeth during one period of shaft rotation. Owing to the similarity between involute gear and involute spline, the random cumulative pitch deviations can be defined as the superposition of sine wave and normal distribution according to the generation mode of tooth indexing errors in the published research about involute spline [32]. Here, the normal distribution is employed to simulate the random characteristic of cumulative pitch deviations, i.e., λ~(µ, σ 2 , where the mean value of cumulative pitch deviations µ = 0 µm, and the standard deviation σ = 3 µm. The schematic diagram of random cumulative pitch deviations and the corresponding pitch deviations are given in Figure 8. The amplitude of random cumulative pitch deviations is set to 15 µm.

Quasi-Static Analysis
For a gear pair with manufacturing errors, the contact patterns will differ under different load conditions. The change of contact patterns will cause the changes of dynamic mesh excitations. In order to further analyze the vibration responses of the HGRS, the quasi-static mesh excitations of the helical gear pair with random cumulative pitch deviations are investigated in this section. LPMS,

Quasi-Static Analysis
For a gear pair with manufacturing errors, the contact patterns will differ under different load conditions. The change of contact patterns will cause the changes of dynamic mesh excitations. In order to further analyze the vibration responses of the HGRS, the quasi-static mesh excitations of the helical gear pair with random cumulative pitch deviations are investigated in this section. LPMS, STE, and CMS of the helical gear pair with random cumulative pitch deviations under different load conditions are given in Figure 9. The load torques are 200 N·m, 350 N·m, 600 N·m, and 1000 N·m, respectively. It can be seen that the cumulative pitch deviations have significant influence on LPMS, STE, and CMS. From Figure 9a, it can be found that the values of LPMS under light load conditions in some mesh cycles are obviously smaller than that under heavy load conditions, and the corresponding curve shapes are notably different from each other. This is because partial contact loss occurs in these mesh cycles under relatively light load conditions. In Figure 9b, the shapes of STE curves under different load conditions are always consistent with each other, but the values become larger in each mesh cycle when the load torque increases. This is because the increased load torque results in larger tooth deflection which contributes to the STE values. From Figure 9c, it can be found that the CMS is related to the load condition. When the load torque is larger than 600 N·m, the curve of CMS remains constant. The changes of LPMS, STE, and CMS can be further explained from the remarkable change of corresponding load distribution. Figure 10 shows the load distributions in the twelfth mesh cycle at different load torques. It can be seen that partial contact loss occurs when the load torque is less than 600 N·m. Moreover, the contact loss region has a tendency to become smaller with the increase of load torque, and the mating gear teeth reaches full contact condition in each mesh cycle when load torque is larger than 600 N·m.      Figure 11 shows the dynamic transmission errors (DTEs) of the helical gear pair with random cumulative pitch deviations under different load conditions. The load torques are also 200 N·m, 350 N·m, 600 N·m, and 1000 N·m, respectively. The input speed is 4000 r/min. It can be seen that the frequency spectrum becomes complex considering random cumulative pitch deviations. In frequency domain of DTEs, the amplitude of shaft frequency is the largest and the sideband frequencies around mesh frequency and its harmonics (MFIHs) appear. When the load torque is smaller than 600 N·m, the amplitudes of MFIHs are less than that of sideband frequencies f m ± i f s (i = 1, 2, . . . , n). With the increase of load torque, the amplitudes of MFIHs will be enhanced. The mesh frequency amplitude become larger than that of sideband frequencies when load torque is higher than 1000 N·m. frequency spectrum becomes complex considering random cumulative pitch deviations. In frequency domain of DTEs, the amplitude of shaft frequency is the largest and the sideband frequencies around mesh frequency and its harmonics (MFIHs) appear. When the load torque is smaller than 600 N·m, the amplitudes of MFIHs are less than that of sideband frequencies ± = ( 1,2,..., ) m s f if i n . With the increase of load torque, the amplitudes of MFIHs will be enhanced. The mesh frequency amplitude become larger than that of sideband frequencies when load torque is higher than 1000 N·m.

Dynamic Responses Analysis
Mesh frequency

Mesh frequency Mesh frequency
Mesh frequency Figure 11. DTEs at different load torques considering random cumulative pitch deviations.
For bearing 1, the vibration acceleration along y direction at different load torques are shown in Figure 12. As similar as DTE, the frequency components of bearing vibration acceleration also become complex. The sideband frequencies ± = ( 1, 2,..., ) m s f if i n appear in the two sides of MFIHs. Some sideband frequencies are larger than mesh frequency when load torque is smaller than 600 N·m. As the load torque increases, mesh frequency amplitude will increase correspondingly. When the load torque is larger than 1000 N·m, the mesh frequency becomes the most predominant. Unlike DTEs, the shaft frequency is not the most predominant frequency in bearing vibration acceleration. For bearing 1, the vibration acceleration along y direction at different load torques are shown in Figure 12. As similar as DTE, the frequency components of bearing vibration acceleration also become complex. The sideband frequencies f m ± i f s (i = 1, 2, . . . , n) appear in the two sides of MFIHs. Some sideband frequencies are larger than mesh frequency when load torque is smaller than 600 N·m. As the load torque increases, mesh frequency amplitude will increase correspondingly. When the load torque is larger than 1000 N·m, the mesh frequency becomes the most predominant. Unlike DTEs, the shaft frequency is not the most predominant frequency in bearing vibration acceleration.

Design of Modification Parameters
It is generally known that the STE fluctuation is usually regarded as an effective standard to evaluate the vibration level of gear pairs. By reducing the material from the tooth top and tooth root, profile modification can change the curve shape and amplitude of transmission errors of a gear pair. When the amplitude of transmission errors is reduced, the vibration of the gear system will become weak. Profile modification has been widely used in the helical gear system [21,[33][34]. Thus, it is necessary to understand the effect of random cumulative pitch deviations on the vibration responses of modified HGRS. Figure 13 shows the design of profile modification for helical gears. Profile modification parameters include modification curve, amount of modification and length of modification. Both driving and driven gear are modified using tip relief. Based on POA, profile modification parameters are defined by the coordinate system XOY, and the potential maximum length of modification is 1 2 B B , max e is the maximum amount of modification, L is the length of modification.

Design of Modification Parameters
It is generally known that the STE fluctuation is usually regarded as an effective standard to evaluate the vibration level of gear pairs. By reducing the material from the tooth top and tooth root, profile modification can change the curve shape and amplitude of transmission errors of a gear pair. When the amplitude of transmission errors is reduced, the vibration of the gear system will become weak. Profile modification has been widely used in the helical gear system [21,33,34]. Thus, it is necessary to understand the effect of random cumulative pitch deviations on the vibration responses of modified HGRS. Figure 13 shows the design of profile modification for helical gears. Profile modification parameters include modification curve, amount of modification and length of modification. Both driving and driven gear are modified using tip relief. Based on POA, profile modification parameters are defined by the coordinate system XOY, and the potential maximum length of modification is B 1 B 2 , e max is the maximum amount of modification, L is the length of modification. Appl. Sci. 2020, 10, x FOR PEER REVIEW 13 of 17 The modification curve adopts quadratic parabola. The equation of profile modification curve can be written as where i x is the distance from the contact point i to the starting position of profile modification, i e is the amount of modification for contact point i.
Because the enhanced LTCA model has very high efficiency, aiming to minimize the STE fluctuation, the optimal profile modification parameters can be obtained by brute global optimization method. This method can clearly show that the change law of STE fluctuation with the modification parameters [17]. The optimization scheme of profile modification is given in Figure 14. It can be observed that there is a narrow region for modification parameters (marked by red line) that reduce STE fluctuation significantly. The STE fluctuation in this region is less than 0.5 μm. From Figure 14, the optimal modification parameters can be obtained very conveniently.  The modification curve adopts quadratic parabola. The equation of profile modification curve can be written as where x i is the distance from the contact point i to the starting position of profile modification, e i is the amount of modification for contact point i.
Because the enhanced LTCA model has very high efficiency, aiming to minimize the STE fluctuation, the optimal profile modification parameters can be obtained by brute global optimization method. This method can clearly show that the change law of STE fluctuation with the modification parameters [17]. The optimization scheme of profile modification is given in Figure 14. It can be observed that there is a narrow region for modification parameters (marked by red line) that reduce STE fluctuation significantly. The STE fluctuation in this region is less than 0.5 µm. From Figure 14, the optimal modification parameters can be obtained very conveniently. The modification curve adopts quadratic parabola. The equation of profile modification curve can be written as where i x is the distance from the contact point i to the starting position of profile modification, i e is the amount of modification for contact point i.
Because the enhanced LTCA model has very high efficiency, aiming to minimize the STE fluctuation, the optimal profile modification parameters can be obtained by brute global optimization method. This method can clearly show that the change law of STE fluctuation with the modification parameters [17]. The optimization scheme of profile modification is given in Figure 14. It can be observed that there is a narrow region for modification parameters (marked by red line) that reduce STE fluctuation significantly. The STE fluctuation in this region is less than 0.5 μm. From Figure 14, the optimal modification parameters can be obtained very conveniently.

Dynamic Responses Analysis
Considering random cumulative pitch deviations, the DTEs of the helical gear pair with and without profile modification are given in Figure 15. For the unmodified gear pair, the amplitude of shaft frequency in DTEs is the largest. Comparing Figure 15a with Figure 15b, it can be found that the amplitudes of MFIHs in DTEs decrease remarkably due to profile modification, but the corresponding sideband frequencies around MFIHs show a slight increasing trend. Considering random cumulative pitch deviations, for bearing 1, the vibration acceleration along y direction are given in Figure 16. As similar as DTEs, the amplitudes of MFIHs decrease significantly, but the amplitudes of sideband frequencies around MFIHs tend to be increased notably.

Dynamic Responses Analysis
Considering random cumulative pitch deviations, the DTEs of the helical gear pair with and without profile modification are given in Figure 15. For the unmodified gear pair, the amplitude of shaft frequency in DTEs is the largest. Comparing Figure 15a with Figure 15b, it can be found that the amplitudes of MFIHs in DTEs decrease remarkably due to profile modification, but the corresponding sideband frequencies around MFIHs show a slight increasing trend. Considering random cumulative pitch deviations, for bearing 1, the vibration acceleration along y direction are given in Figure 16. As similar as DTEs, the amplitudes of MFIHs decrease significantly, but the amplitudes of sideband frequencies around MFIHs tend to be increased notably. Considering different amplitudes of random cumulative pitch deviations, the root mean square (RMS) values of DTEs of the helical gear pair with and without profile modification are shown in Figure 17. The input speed changes from 25 r/min to 8000 r/min with a step of 25 r/min. For case 1 to case 4, the amplitudes of random cumulative pitch deviations are 0 μm, 5 μm, 10 μm, and 15 μm, respectively. From Figure 17a, it can be observed that the RMS values of DTEs of modified helical gear pair decrease significantly in all input speed conditions when compared with the unmodified gear pair. From Figure 17b, it can be found that profile modification shows a well effectiveness in reducing vibration of the HGRS when the input speed is less than about 5500 rpm, but the RMS values of DTEs of the modified gear pair are larger than that of unmodified one in some relatively higher input speed conditions. This is because the effects of random cumulative pitch deviations on dynamic

Dynamic Responses Analysis
Considering random cumulative pitch deviations, the DTEs of the helical gear pair with and without profile modification are given in Figure 15. For the unmodified gear pair, the amplitude of shaft frequency in DTEs is the largest. Comparing Figure 15a with Figure 15b, it can be found that the amplitudes of MFIHs in DTEs decrease remarkably due to profile modification, but the corresponding sideband frequencies around MFIHs show a slight increasing trend. Considering random cumulative pitch deviations, for bearing 1, the vibration acceleration along y direction are given in Figure 16. As similar as DTEs, the amplitudes of MFIHs decrease significantly, but the amplitudes of sideband frequencies around MFIHs tend to be increased notably. Considering different amplitudes of random cumulative pitch deviations, the root mean square (RMS) values of DTEs of the helical gear pair with and without profile modification are shown in Figure 17. The input speed changes from 25 r/min to 8000 r/min with a step of 25 r/min. For case 1 to case 4, the amplitudes of random cumulative pitch deviations are 0 μm, 5 μm, 10 μm, and 15 μm, respectively. From Figure 17a, it can be observed that the RMS values of DTEs of modified helical gear pair decrease significantly in all input speed conditions when compared with the unmodified gear pair. From Figure 17b, it can be found that profile modification shows a well effectiveness in reducing vibration of the HGRS when the input speed is less than about 5500 rpm, but the RMS values of DTEs of the modified gear pair are larger than that of unmodified one in some relatively higher input speed conditions. This is because the effects of random cumulative pitch deviations on dynamic Considering different amplitudes of random cumulative pitch deviations, the root mean square (RMS) values of DTEs of the helical gear pair with and without profile modification are shown in Figure 17. The input speed changes from 25 r/min to 8000 r/min with a step of 25 r/min. For case 1 to case 4, the amplitudes of random cumulative pitch deviations are 0 µm, 5 µm, 10 µm, and 15 µm, respectively. From Figure 17a, it can be observed that the RMS values of DTEs of modified helical gear pair decrease significantly in all input speed conditions when compared with the unmodified gear pair. From Figure 17b, it can be found that profile modification shows a well effectiveness in reducing vibration of the HGRS when the input speed is less than about 5500 rpm, but the RMS values of DTEs of the modified gear pair are larger than that of unmodified one in some relatively higher input speed conditions. This is because the effects of random cumulative pitch deviations on dynamic responses are not the same in different input speeds due to natural frequencies of gear systems and frequency components of mesh excitations. As the amplitude of random cumulative pitch deviations increases, the influences of cumulative pitch deviations on the effectiveness of profile modification in reducing system vibration become more obvious. When the amplitude of random cumulative pitch deviations is larger than 10 µm, the vibration of the modified HGRS becomes stronger than that of the unmodified one in most input speed conditions, as shown in Figure 17c,d. Thus, as a long-period component of gear errors, cumulative pitch deviations not only affect the kinematic accuracy of gear pairs, but also have significant effect on vibration responses of HGRS.
Appl. Sci. 2020, 10, x FOR PEER REVIEW 15 of 17 responses are not the same in different input speeds due to natural frequencies of gear systems and frequency components of mesh excitations. As the amplitude of random cumulative pitch deviations increases, the influences of cumulative pitch deviations on the effectiveness of profile modification in reducing system vibration become more obvious. When the amplitude of random cumulative pitch deviations is larger than 10 μm, the vibration of the modified HGRS becomes stronger than that of the unmodified one in most input speed conditions, as shown in Figure 17c,d. Thus, as a long-period component of gear errors, cumulative pitch deviations not only affect the kinematic accuracy of gear pairs, but also have significant effect on vibration responses of HGRS.

Conclusions
In this paper, the influences of random cumulative pitch deviations on mesh excitations and vibration responses of HGRS are studied. The LPMS, STE, and CMS of helical gear pairs are determined using an enhanced LTCA model. A dynamic model with multi DOF is employed to predict the vibration responses of HGRS. The effect of random cumulative pitch deviations on vibration responses of modified HGRS in time and frequency domain are also investigated in detail. The main findings are summarized as follows: (1) considering random cumulative pitch deviations, the curve shapes of STE and CMS become irregular. When the load torque is relatively low, the LPMS in some mesh cycles will decrease due to partial contact loss. As the load torque increases, the actual contact region will become larger. When the contact patterns in each mesh cycle extend to full tooth surface, the curve shape of LPMS will become regular; (2) vibration modulation phenomenon can be observed in dynamic responses of HGRS with random cumulative pitch deviation, and the random characteristic of cumulative pitch deviations makes the vibration frequency spectrum complex. When HGRS is under light load conditions, the amplitudes of some sideband frequencies are larger than that of MFIHs because

Conclusions
In this paper, the influences of random cumulative pitch deviations on mesh excitations and vibration responses of HGRS are studied. The LPMS, STE, and CMS of helical gear pairs are determined using an enhanced LTCA model. A dynamic model with multi DOF is employed to predict the vibration responses of HGRS. The effect of random cumulative pitch deviations on vibration responses of modified HGRS in time and frequency domain are also investigated in detail. The main findings are summarized as follows: (1) considering random cumulative pitch deviations, the curve shapes of STE and CMS become irregular. When the load torque is relatively low, the LPMS in some mesh cycles will decrease due to partial contact loss. As the load torque increases, the actual contact region will become larger. When the contact patterns in each mesh cycle extend to full tooth surface, the curve shape of LPMS will become regular; (2) vibration modulation phenomenon can be observed in dynamic responses of HGRS with random cumulative pitch deviation, and the random characteristic of cumulative pitch deviations makes the vibration frequency spectrum complex. When HGRS is under light load conditions, the amplitudes of some sideband frequencies are larger than that of MFIHs because of relatively high contact ratio. With the increase of load torque, the amplitudes of MFIHs will increase and become larger than that of sideband frequencies; (3) profile modification can reduce the amplitudes of MFIHs in DTEs and bear vibration acceleration.
However, with the increase of the amplitude of cumulative pitch deviations, the effectiveness of profile modification in reducing system vibration will become weaker. In some working conditions, profile modification will lead to stronger system vibration for HGRS, owing to the increased amplitudes of sideband frequencies.