Relationship between the engaging force of planetary gear train and the position-correlated modal properties

: Planetary gear train (PGT) is widely used in a variety of fields such as helicopters and aircraft engines. The failure of meshing gear often happens in the PGT due to the working condition of high speed, variable load, high-frequency excitation, assembly error etc. However, the causes of gear failure are still insufficiency, especially in the condition of high speed and high-frequency excitation. In this work, the engaging force between gear meshing pair of the PGT under high-frequency excitation is studied using the proposed position-correlated modal properties calculation method, which is established by incorporating the effect of meshing position and meshing phase difference of each contact pair into the free vibration model to study the modal properties of the PGT at different meshing positions. The corresponding engaging force responses based on the position-correlated modal properties are investigated. The simulation results show that the higher-order natural frequencies of the system are greatly affected by the meshing position. The peaks of engaging force occur at meshing positions where the natural frequencies are equal to the excitation frequency, which can be the potential cause of the gear damage.


Introduction
Planetary gear train (PGT) is widely used in mechanical transmission due to its higher transmission ratio and highly compact structure. The engaging force and modal properties of the PGT are the fundamental factors determining the performance of the mechanical transmission system and the machine structure.
The engaging force of the PGT is complex due to its multiple degrees of freedom (DOFs) and meshing phase difference between the multiple meshing gear pair. On the other hand, the engaging force is affected by the natural frequencies and vibration modes of the planetary gear system since the engaging force will be greatly increased when the system resonant is excited. Parker and Ambarisha [1,2] used Fourier series method to study the engaging force of the PGT, which is further used to predict whether a particular type of mode is excited.
Due to the significance of modal properties of the PGT, considerable attention has been focused on the free vibration model of the PGT and the system modal properties. A purely torsional free vibration models of the PGT and compound PGT were established by Kahraman [3,4] to analyse its modal properties. Lin and Parker [5] studied the modal properties of the torsionaltranslational model and summarised three types of free vibration modes as rotational mode, translational mode, and planet mode. Based on compound multi-stage planetary gear sets, Kiracofe and Parker [6] extended the analytical models to study the vibration characteristics of the system. Guo and Parker [7] also used a purely rotational model to analyse the modal properties of compound planetary gear system. The structured modal properties of the helical PGT are categorised in [8]. Sensitivity analysis of natural frequencies and vibration modes of planetary gear system is derived by Lin and Guo in Ref [9,10]. Wu and Parker [11] determined modal properties of planetary gears with an elastic continuum ring gear. A coupled concentrated parameter model of a multi-stage PGT is established by Sun [12] to analyse the coupled modes in multi-stage PGT. Ambarisha and Parker [2] used modal force and mesh phasing to study the suppression phenomenon of planet mode response in planetary gear dynamics. Natural frequency veers and clusters in planetary gear vibration are observed in [13,14]. Zhang [15] studied the natural frequencies and mode shapes of a two-stage closed-form PGT, and parameters' influence on modal characteristics is analysed. Kim [16] partitioned the boundary conditions of the transfer matrix into known and unknown values to generate a concentrated transfer matrix and a latent equation, solved the eigenvalue problem in the lambda matrix, and studied the mode shape type of the PGT.
Investigations of modal properties of the PGT have been carried out extensively, but litter work is done on the effect of the meshing position and meshing phase difference of each contact pair on the system modal properties and the further effect on engaging force. In practical meshing state of planetary gears, the meshing stiffness is unique at different contact points on the planet meshing line, which means that the system has unique property at different meshing positions. Results in [9,13] show that the variation of meshing stiffness has a remarkable influence on the higher-order natural frequency of the system. With this property, if excitation frequency happens to close to the natural frequencies in a specific position, the amplitude of meshing force could significantly increase, which can further increase the probability of tooth damage. Therefore, the study of the modal properties based on the meshing position is necessary to the further study of the engaging force response.
On the other hand, the PGT, which is most commonly used in product applications, is the one with equally spaced planets and sequentially phased gear meshes since it was shown to result in lower vibration and noise levels. In this case, the meshing phase difference exists in each contact gear and the engaging status of each gear is different. Lin and Parker [9] mentioned that different meshing stiffness between planets can be caused by the factors of different number of teeth in contact, manufacturing variation and assembly errors, which result in the mistuning of the engaging status of each planet mesh and affect the modal properties of the system. Due to the meshing phase difference between planets, the meshing stiffness of each planet engaging with the central gear is different and the phenomenon of mistuning meshing status will always exist, which also influences the modal properties. The influence of meshing phase on the PGT modal property should be also considered for the further study of engaging force response.
Due to the above reasons, in order to study the engaging force response of the system, it is essential to establish a free vibration model, which takes the meshing position and meshing phase of each planet into consideration to comprehensively study the effect of meshing position and meshing phase on the modal properties of the PGT. Based on the position-correlated modal properties, a new insight for the analysis of the engaging force response can be put forward.
In this paper, a purely torsional free vibration model considered the meshing position is established. The influence of meshing position on the modal properties is studied. The meshing phase difference is also taken into account to investigate the distribution regulation of the natural frequencies in a position interval to lay a foundation for the analysis of engaging force under high-frequency excitation. The results show that the modal properties of the system show disparate characteristics at different meshing positions and the system has abundant modal properties under higher-order natural frequencies. The number and interval of resonant positions in one meshing cycle depend on the meshing phase. Due to the relatively wide distribution region of higher-order frequencies in a meshing position interval, the regions of excitation frequency which can arouse high-amplitude engaging force are wider in high frequencies. These results provide another insight into the engaging force analysis of the PGT under high-frequency excitation.
The rest of the paper is structured as follows. Section 2 is the free vibration model description of the PGT. Section 3 is the derivation of the relation between modal properties of the different meshing positions and meshing phase difference. Section 4 is the discussion on the simulation results obtained and analysis of engaging force response which is affected by the positioncorrelated modal properties. Concluding remarks are drawn in Section 5.

System free vibration model
To analyse the modal properties of the system incorporating the meshing position and meshing phase, only rotational DOFs are taken into consideration and translational motions of each component are eliminated in the model. The lumped-parameter model of the planetary gear system is shown in Fig. 1. The ring gear is held stationary, and the mass and moment of inertia of each planet gear are identical. The planet gear rotates around the sun gear and ring gear with mesh stiffness k spn and k rpn . k su , k ru , and k cu are the torsional stiffness of sun, ring gear, and carrier, which are constant. Sun, planet, ring gear, and carrier vibrate in rotational direction with deformations of u s , u n , u c , u r , and the counterclockwise circumferential positioning angles of nth planet relative to first planet are denoted by ψ n = 2π(n − 1)/N, ψ 1 = 0.
With the change of circumferential position of sun gear, the meshing state of the system is also changed, without any loss of generality; θ is defined as the circumferential position of the sun gear with respect to initial position to describe the mesh status of each planet; and the initial position starts from horizontal direction, which is shown in Fig. 2a, θ 1 = 0. The circumferential position of sun gear varying from position θ 1 to position θ 2 is illustrated in Fig. 2b, since the meshing status in θ 2 is different from the meshing status in θ 1 , and the modal properties of these two positions are distinct.
Introduce parameter κ xp1 θ x = s, r to describe the mesh stiffness variation function of the x/n th planet mesh at the position θ of sun gear, i.e.
In (1), γ sr denotes the relative meshing phase difference between the mesh, Θ m denotes the circumferential angle interval corresponding to one complete engagement cycle which yields where z s denotes the teeth of sun gear, and z r denotes the teeth of the ring gear. The mesh stiffness of the first sun/planet k sp1 θ and nth sun/ planet gear k spn θ is identical and consequently ring-planet mesh stiffness due to the symmetry distribution of the planet gears and the same parameters of all planet gears. The only difference is that there is a meshing phase difference between k sp1 θ and k spn θ (also k r p1 θ and k r pn θ ) because the mesh status of each planet gear is not in phase with each other. The expressions of k spn θ and k r pn θ can be written as follows: where γ xpn x = s, r denotes the meshing phase between the x/nth planet mesh and the x/1st planet mesh, which can be defined as the following equation: (see (4)) . The assembly condition of equally spaced planetary gear system has Substituting (5) into (4), the expression of γ xpn x = s, r can be written as (see (6)) . Therefore, for equally spaced gear system, the relationship between γ spn and γ r pn satisfies At the circumferential position of sun gear θ, The undamped free vibration equations of the planetary gear system can be written as follows: (see (8)) where N is the number of planetary gear. = u s + u n − u c cos α s , δ r pn = u r − u n − u c cos α r , where α s and α r are pressure angles of sun gear/planet gear mesh and ring gear/ planet gear mesh.
Equation (8) can be rewritten in the form of matrices as where M is the inertia matrix, K m θ and K b are stiffness matrices, u is rotational motion matrix, and they are expressed as . K θm [θ should be formatted as superscript] [m should be formatted as subscript] is shown in Fig. 3. Assuming the motion of all the components in the planetary gear system under free vibration is harmonic vibration, the associated eigenvalue for (9) can be written as where ω i is the ith order of the natural frequency at the position of sun gear θ, ϕ i denotes the corresponding vibration mode under ω i with the form： By solving the characteristic equation (10), the modal properties of a given position θ can be determined.

Effect of meshing phase on modal properties
Due to the existence of the meshing phase difference between two adjacent planet gear meshes, the meshing state between two planet gears has a phase lag or lead, and the meshing phase difference γ (lag or lead) of two adjacent planet meshes is identical. Therefore, the meshing state of the nth planet mesh at the position of θ is coincident with the meshing state of the (n + 1)th planet at position of θ + γΘ m , which means that the system modal properties at the position θ and θ + γΘ m are correlated. The mathematical derivation process is shown as follows.
With counterclockwise meshed planet, based on the expression of meshing phase shown in (6), the (n + 1)th planet teeth mesh lags n th planet mesh by a phase angle γ′ which is defined by the following equation: where the subscription x denotes s(sun gear) and r(ring gear), n stands for the serial number of nth planet, N denotes the number of the planet gears, mod(a, 1) stands for the remainder after division of a by 1.
Introduce the parameter of phase difference γ between two adjacent planets' tooth meshing, which is expressed as For 0 < γ′ ≤ 0.5, the (n + 1)th planet tooth mesh lags nth planet mesh by a phase γΘ m . For 0.5 < γ′ < 1, the (n + 1)th planet tooth mesh advances nth planet mesh by a phase γΘ m . According to the identical meshing phase difference between two adjacent planetary gears, the relation of mesh stiffness of two planet gears at θ and θ + n′γΘ m satisfies the following: where n′ denotes an arbitrary positive integer.
In (14), if n + n′ > N, it has where floor(x) denotes round x to the nearest integer which is ≤ x, and ceil(x) denotes round x to the nearest integer which is ≥ x.
Combining (14) and (16), (see (17)) . Substituting (17) into (3), the meshing stiffness can be expressed as follows: To distinctly describe the variation elements in stiffness matrix, k θ (R, C) is introduced to denote the element in K m θ at the Rth row and Cth column.
For R ∈ (1, 2, 3) and C ∈ (1, 2, 3), with (17) and constant pressure angles α s and α r , the relationship of elements can be expressed as follows: Equation (19) indicates that the 3 × 3 sub-matrix at the upper-left corner of the matrix K m θ and K m θ + n′γΘ m is identical, k R, C periodically changes 1/γ times over one meshing period.
For R ∉ (1, 2, 3)andC ∉ (1, 2, 3), contrast θ to θ + n′γΘ m , the position of matrix elements shifts n′ units length in a certain direction. The mathematical expression of this process is shown in (20)(21)(22).  (22)) The k θ R, C and k θ + n′γΘ m R, C determined by (19)- (22) are substituted into (10) to solve the characteristic eigenvalues at the position θ and θ + n′γΘ m . It is noted that the coordinate transformation of the elements in K m is diagonally symmetric, and the eigenvalues k xp1 θ + k xp2 θ + ⋯ + k xp N − 1 θ + k xpN θ = k xp2 θ + γΘ m + k xp3 θ + γΘ m ⋯ + k xpN θ + γΘ m + k xp1 θ + γΘ m = k xp3 θ + 2γΘ m + k xp4 θ + 2γΘ m ⋯ + k xp1 θ + 2γΘ m + k xp2 θ + 2γΘ m = ⋯ 0 < γ′ < 0.5 There are other n positions with the same natural frequency as the θ in this interval θ 0 , θ 0 + Θ m : According to eigenvectors analysis, the eigenvectors of first three elements of each order (corresponding to the free vibration modes of sun gear, carrier, and ring gear) remain unchanged while the other elements (corresponding to the free vibration modes of planet gears) shift downward 0 < γ′ ≤ 0.5 or upward 0.5 < γ′ < 1 by a length of n′ units. There is a relation of the free vibration modes between two specific meshing positions. The free vibration modes of sun gear, ring gear, and carrier at the positions θ and θ + n′γΘ m are coincident with each other. Furthermore, the vibration mode of nth planet at θ′ is consistent with the free vibration mode of (n + 1)th planet 0 < γ′ ≤ 0.5 or (n-1)th planet 0.5 < γ′ < 1 at θ + n′γΘ m . Based on the relation between the meshing phase and modal properties, an important implication is that there are 1/γ positions in one meshing cycle with the identical natural frequencies, which means that if the excitation frequency happens to near natural frequency of a one position, there may appear 1/γ peaks of engaging force in one meshing cycle.

Simulation and discussions
The planetary gear system with fixed ring gear and counterclockwise meshed planet is taken as the example. The parameters are shown in Table 1.
According to (6), the meshing phases of sun/planet gear and ring/planet gear are as follows:   The (n + 1)th planet teeth mesh lags nth planet teeth mesh by a phase angle γ′ = 1/4, which indicates that the meshing phase difference of two adjacent planetary gears is equal to 1/4 γ = 1/4 . The relative phase γ rs is determined as γ rs = 0 using the method in [17].
To obtain the position-correlated modal properties of the system, the position-correlated meshing stiffness is calculated firstly. Supposing the material of contacting teeth is isotropic, the single tooth meshing stiffness can be expressed as follows: (see (26)) where k b , k s , and k a denote the bending stiffness, shear stiffness, and axial compressive stiffness of the single tooth, respectively. k h denotes Hertz contact stiffness, and k f 2 is the stiffness with consideration of gear fillet-foundation deflection. This stiffness is calculated using the potential energy method based on cantilever beam theory [18][19][20][21][22][23]. The meshing stiffness of sunplanet gear mesh and ring-planet gear mesh can be obtained. Taking meshing phase between planets into account, with θ 0 representing the position that one pair of meshing teeth between sun gear and first planet is just going into meshing, the correspondence mesh stiffness k xpn of internal and external engagement of each planet gear at different positions over one meshing cycle is illustrated in Fig. 4. It can be noticed that the (n + 1)th planet tooth mesh lags nth planet mesh by a phase Θ m /4 and with the existence of meshing phase difference, and the meshing stiffness of each planet is different at any given positions.

Analysis of position-correlated natural frequencies
By solving (10) in critical positions, the distribution of positioncorrelated natural frequencies in interval θ 0 , θ 0 + Θ m can be obtained and the modal properties of the system in these mesh positions can be determined, as shown in Fig. 5. The positionindependent natural frequencies determined by ignoring the meshing position of meshing gears with meshing stiffness equalling to their mean values, k sp1 = k sp2 = ⋯ = k spN = k sp mean and k r p1 = k r p2 = ⋯ = k r pN = k r p mean , are also illustrated in Fig. 5.
The minimum natural frequencies at any positions are zero corresponding to the rigid body motion of free vibration mode, due to the positive semi-definite stiffness matrix K m θ . For higher-order frequencies calculated by the proposed method, as shown in Fig. 5, there is a dash-dotted line region, called 'natural frequency concentrated region' where the natural frequencies of third to sixth order distribute very intensively and densely and the position-independent natural frequencies of third to sixth happen to fall in this region.
In the operation process of the planetary gear system, the meshing state of each planet is changing periodically. If the excitation frequency happens to approach the position-correlated frequency of a certain mesh position, high-amplitude meshing force would be aroused with high probability at this position. So excitation frequency should better be kept away from this 'natural frequency concentrated region' to avoid excessive impact.
The position-correlated natural frequencies and positionindependent natural frequencies of each order in position interval θ 0 , θ 0 + Θ m are separately shown in Fig. 6. For position-correlated natural frequencies, it can be noticed that the modal of the system shows different characteristics at different meshing positions and the natural frequencies regularly distribute in this position interval. Four identical natural frequencies can be observed in the studied position interval. The number is equal to the theoretical number calculated using (12), (13), (24): The interval between two adjacent position points is the same which is equal to Θ m γ = Θ m /4: A significant implication is that if the excitation frequency just falls in this 'natural frequency concentrated region' 1/γ, peaks of engaging force could be noticed and the system will experience 1/γ times periodic impact due to the meshing phase. Fig. 6 also illustrates that the fluctuation range for the positioncorrelated frequency of the second order is limited in a small range (<7 Hz) and the value is almost the same as the position of independent natural frequency. The frequency variation range of sixth-order position-correlated frequency is also small (<20 Hz), but there is a value deviation between the two methods. For other orders (third, fourth, fifth, seventh), the system natural frequencies are dramatically affected by the location of the mesh position, and there is a very large variation range (>500 Hz) for the positioncorrelated natural frequencies.
A distinct loci change phenomenon of the natural frequencies is obviously shown in Fig. 6. The position-correlated natural frequencies ω n n = 2, …, N vary abruptly near the position θ = a, b, c etc. These positions are in accordance with the saltation points of the meshing stiffness, as the points shown in Fig. 4a corresponds to the termination position of the double teeth mesh area of third planet/sun gear mesh, b corresponds to the termination position of the double-teeth mesh area of second planet/sun gear mesh, c corresponds to the termination position of the single tooth mesh area of second planet/sun gear mesh and second planet/ring gear mesh. This is due to the changed pair of contact teeth, which leads to a dramatic change of mesh stiffness and the further large variation of modal properties. This phenomenon is especially obvious in higher-order natural frequency (>7000 Hz), which suggests that the differing number of teeth in contact has great influence on higher-order natural frequencies.
In order to demonstrate the effect of meshing position on the system modal properties more generally, another planetary gear system with torsional and translational DOFs for each component is further studied. The parameters in Table 1 are also employed for the calculation. Fig. 7 illustrates the natural frequencies versus the position of the sun gear in interval θ 0 , θ 0 + Θ m with the translational freedom included using both position-correlated mesh stiffness and mean value mesh stiffness. After considering the translational DOFs, the distribution of natural frequencies is denser than the results calculated by rotational model, but the location intervals are broadly consistent with each other. As can be observed in Fig. 7, the sudden change property of mesh stiffness can largely cause the value oscillation of the higher-order natural frequencies. Similar to that presented in Figs. 6 and 7, the two methods predict almost the same natural frequencies for low-order frequencies. However, a difference is observed between the two methods for higher-order natural frequencies. The higher-order position-correlated natural frequencies are highly dependent on the meshing position; especially, the differing number of teeth in contact of each planet has large influence on the natural frequencies.

Analysis of engaging force based on position-correlated modal properties
As the natural frequencies of the system are distinct at different meshing positions, the following situation is possible: the system excitation frequency is different from the natural frequency at one mesh position θ 1 ′, and the excessive impact will not occur. However, in the real operation process of the PGT, because of the variation of the meshing position, in position θ 2 ′, the corresponding changed natural frequency approaches to the excitation frequency ω m , the amplitude of the engaging force will rapidly increase, and the excessive impact occurs which will result in the increase of the root mean square (RMS) value of engaging force. The frequent impact will further increase the probability of tooth damage.
In order to vindicate the assumption, engaging force response of the example spur planetary gear set (Tab.1) is simulated under the excitation frequency 8700 Hz with sun gear input (3000 rpm) and carrier output. The engaging force responses of the planet/ring gear mesh in the position interval [30θ 0 , 30θ 0 + Θ m ] are shown in Fig. 8a. It can be noticed that the peaks of engaging force occur at positions where the natural frequencies are equal to the excitation frequency 8700 Hz, as it is shown in Fig. 8b. In addition, the number of these positions is equal to 1/γ and the distance between two adjacent positions is equal to γΘ m , which indicates that the system will experience 1/γ times impact due to the meshing phase. These simulated results agree quite well with the theoretical derivation in Chapter 3, which indicates that the position-correlated modal properties analysis can provide a useful guidance for the analysis of engaging force response under high-frequency excitation. Fig. 9 shows the dynamic response of the RMS value of ring/ first planet engaging force F r/ p (rms) under the varying mesh frequency ω m . An obvious peak of F r/ p (rms) can be observed at the mesh frequency ω m = 2017 Hz. In the region around ω m ≃ 2107 Hz, F r/ p (rms) rises rapidly firstly, as soon as it reaches the resonance peak, and F r/ p (rms) falls off sharply within a short time. It can be noted that the excitation frequency, which corresponds to simulation peak point, is in good agreement with the two kinds of natural frequency analysed incorporating position-independent mesh stiffness (ω n = 2017.7 Hz) or position-correlated mesh stiffness (ω n ∈ 2013 Hz, 2010 Hz ). For meshing frequency, ω m moves to high-frequency region of 7950-9410 Hz, and the corresponding F r/ p (rms) is significantly larger than other regions. The excitation frequency ω m in this region can be regarded as the frequency that can excite excessive impact of the system. This excessive impact region is in accordance with the natural frequency concentrated region shown in Fig. 5. The results show that the natural frequency concentrated region is the potential cause of gear tooth damage of the PGT under the high-frequency excitation. If the excitation frequency or its multiple frequency just falls in this range, the amplitude of engaging force is excessively increased, which may result in the damage of the teeth. According to the results obtained, to avoid the excessive meshing force, the excitation frequency ω m and its lth harmonics lω m should better avoid falling in this natural frequency concentrated region. These results provide a useful insight into the fault diagnosis of the PGT under high-frequency excitation.

Conclusions
This work analytically investigates the modal properties of torsional vibration model for planetary gear transmission based on meshing position. Effects of position-correlated mesh stiffness and meshing phase between two adjacent planets on the modal properties of the system are presented. Depending on the positioncorrelated modal properties, the corresponding engaging force response with high-frequency excitation is investigated. The main findings are as follows: (i) The meshing position has influence on the modal properties, and the system modal shows disparate characteristics at different meshing positions, especially the higher-order natural frequencies of the system. In the case of analysing the engaging force response under the high excitation frequencies, the effect of meshing position should be taken into consideration. (ii) There is a relation between the position-dependent modal properties and the engaging force response. The peaks of engaging force occur at meshing positions where the natural frequencies are equal to the excitation frequency. The number and interval of these positions in one meshing cycle depend on the meshing phase. (iii) There exists a 'natural frequency concentrated region' in higher-order frequency range, which is the potential cause of gear tooth damage of the PGT under the high-frequency excitation. If the excitation frequency or its multiple frequency just falls in this range, the amplitude of engaging force is excessively increased, and the excitation frequency and its harmonics should avoid falling in this region to prevent the excessive engaging force.