Planetary gearbox with localised bearings and gears faults: simulation and time/frequency analysis

Planets bearings of planetary gear sets exhibit high rate of failure; detection of these faults which may result in catastrophic breakdowns have always been challenging. The objective of this paper is to investigate the planetary gears vibration properties in healthy and faulty conditions. To seek this goal a previously proposed lumped parameter model (LPM) of planetary gear trains is integrated with a more comprehensive bearing model. This modified LPM includes time varying gear mesh and bearing stiffness and also nonlinear bearing stiffness due to the assumption of Hertzian contact between the rollers/balls and races. The proposed model is completely general and accepts any inner/outer race bearing defect location and profile in addition to its original capacity of modelling cracks and spalls of gears; therefore, various combinations of gears and bearing defects are also applicable. The model is exploited to attain the dynamic response of the system in order to identify and analyze localized faults signatures for inner and outer races as well as rolling elements of planets bearings. Moreover, bearing defect frequencies of inner/outer race and ball/roller and also their sidebands are discussed thoroughly. Finally, frequency response of the system for different sizes of planets bearing faults are compared and statistical diagnostic algorithms are tested to investigate faults presence and growth.


Introduction
Planetary gears, also recognized as epicyclic gears, are extensively used power transmission elements in numerous fields such as automotive, aerospace, wind turbines and marine applications. They have several benefits including compactness, high torque to weight ratio, high efficiency, multiple gear ratios and reduced noise in comparison with fixed-shaft gearboxes. Therefore, investigating planetary gear noise and vibration in healthy and faulty conditions is crucial to keep them functional and also to avoid any machinery breakdown as a result of a partial failure.
Mathematical modeling is an advantageous approach to scrutinize characteristics of mechanical systems. It gives a good understanding of structure dynamic characteristics; it is reasonably accurate and suitable for evaluations during design stages. In this regard mathematical models such as lumped parameter models (LPMs) have been vastly used to study modal properties and also vibration signals of planetary gear trains.
Cunliffe et al. [1] developed a two dimensional mathematical model for a planetary gearbox with thirteen degrees of freedom and studied the natural frequencies and vibration modes of the system. Botman [2] analyzed the effect of planet-bearing stiffness and rotation of the carrier on the natural frequencies of in-plane vibration of a single stage spur planetary gearbox with eighteen degrees of freedom. Frater et al. [3] extended Botman's model to take into account the unequal mesh stiffness between the sunplanet and ring-planet meshes and also stiffness alternation due to variation of number of teeth in contact. Kahraman [4] developed a nonlinear timevarying dynamic model of a planetary gear set which includes tooth separations and mesh stiffness fluctuations. He subsequently derived closed form expressions for natural modes of a planetary gear by using a purely torsional model [5]. Lin and Parker [6] developed an analytical model and rigorously investigated main properties of natural frequencies and vibration modes of a general planetary gear system with equally spaced planets which is applicable to various gears configurations. They later studied the free vibration of this set with unequally spaced planets [7]. Sun and Hu [8] developed a lateral-torsional coupled nonlinear dynamic model for a planetary gear system with multiple clearances. Ambarisha and Parker [9] used both the lumped parameter and finite element models to examine the tooth separations nonlinear dynamic behavior of spur planetary gears. Guo and Parker [10] introduced tooth wedging, tooth contact loss and bearing clearance into a lumped parameter model and investigated the interplay between tooth wedging and bearing clearance. They also [11] discussed nonlinear behavior, bifurcations, and chaos caused by bearing clearance as well as interaction between tooth separation and bearing clearance.
Simulation has also been used to model gear faults due to its inherent advantages such as deeper insight into the signatures of different faults on the system vibration signals, the reduction of number of experiments and the testing of fault diagnostic and prognostic techniques. Chaari et al. [12] studied the influence of tooth defects (pitting and crack) on the response of the system by modifying the gear mesh stiffness. Chen and Shao [13][14][15] presented a series of papers on the effects of gears cracks on planetary gear dynamics. They investigated rigid ring gear tooth crack [13], sun and planet gears teeth cracks with different sizes and inclination angles [14]. Moreover, Chen et al. [15] studied the influence of flexible ring gear rim and cracked tooth on mesh stiffness and dynamic features of a planetary gearbox. Wu and Parker [16] derived the characteristic modal properties of planetary gears having equally spaced planets and an elastic continuum ring gear. Zhang et al. [17] calculated natural frequencies and corresponding vibration modes of a planetary gearbox with flexible ring by dividing the continuum ring into finite rigid sectors. To determine the vibration signatures of localized planet-bearing faults Jian [18] developed an analytical model including flexible ring gear. This paper develops a lumped parameter model to investigate the gears and bearings interaction of a planetary gear train in presence of faults. A 18-DOFs model of planetary gearbox, based on [19], is combined with a comprehensive bearing model. The resulting set of equations has the capability of simulating the behaviour of the system for different sizes, locations and profiles of defects, both on gears and on bearings, and also includes non-linear effects due to the Hertzian contact assumption. The purpose is to build a numerical tool where faults, different in type, size and position, can be implemented to produce time domain signals. The aim is not to exactly reproduce the behaviour of a real gearbox, or even of a test rig, but to produce numerical simulations of faults which are and important step in testing and developing diagnostic techniques. Numerically generated signals, with specific characteristics directly related to particular defects, can provide an efficient and very economical way to verify and enhance the capabilities of data processing methods. In real life, the signature of the defect can in fact be masked by external noise, making it impossible to figure out what the real contribution of the defect on the measured output is; and even more difficult is to have at disposal a test rig with known and controlled defects. Furthermore, the paper comprehensively discusses the frequency components of signals associated to planet-bearing defects for inner/outer race and rolling element and sources of sidebands around bearing damage frequencies. Finally, frequency analysis and statistical features are implemented on The primary aim of this simulation is to calculate vibration signals of each gear of a planetary gearbox in presence of defects on bearings and gears. The twodimensional lumped parameter models by Lin and Parker [6] and Liang et al. [19] are the basis for the present discussion. There are some distinctions between these two models: first, configurations of planets coordinates are different; second, damping is not considered in the Lin and Parker [6] model. Whether planets deflections are described cartesian coordinates [19] or in polar coordinates [6], the dynamic behavior of this gearbox will not change. The coordinates in this paper are based on Ref. [19] which simplify the final expression of the equations of motion.
A 2D lumped parameter model for a spur planetary gearbox is illustrated in Fig. 1. Each gear (sun, ring and planets) and carrier are considered as rigid bodies with three degrees of freedom (DOFs)-one rotational (h) and two transverse motions in the x and y directions. The resulting LPM has 3ðN þ 3Þ DOFs where N is number of planets. The degrees of freedom of this system are: [x c ; y c ; h c ; x r ; y r ; h r ; x s ; y s ; h s ; x pj ; y pj ; h pj ] where c, r, s and pj; j ¼ 1; : : :; N are assigned to the carrier, ring, sun and planet gears, respectively. The flexible gear teeth contacts are modelled by springs and dampers acting along the gear line of action and tooth contact loss is assumed not to occur. The gears translational displacements are calculated with respect to a rotating frame of reference Oxy with origin O, the center of the planetary gear set. This frame is attached to the carrier and rotates with the same angular speed as carrier. The angular rotation h is defined in OXY reference system which is fixed and is not rotating. Furthermore, all gears are assumed perfect without manufacturing and mounting errors. Gears and carriers are also considered free of eccentricities and roundness errors. According to Jian [18], the flexibility of the ring gear influences the relative amplitudes of the sidebands in a fault signature and higher-order sidebands disappear when ring thickness increases. The ring gear herewith described is not deformable but a development of the model is planned so to connect the ring gear, the bearing of the sun gear and the bearing of the carrier to a flexible gearbox casing, to simulate an elastic support and the path from the defect to the measurement point.

Time-varying mesh stiffness
The condition in which gear teeth are in contact varies as they rotate. For a contact ratio lower than two, the number of teeth pairs in contact periodically changes from one to two and this causes a time variation of mesh stiffness, which is the main source of vibration in gearboxes. Mesh stiffness also changes with the contact positions of gear teeth. Teeth pairs enter and exit the mesh constantly and therefore the assumption of equivalent characteristics for every tooth leads to periodic time varying mesh stiffness. Yang and Lin [20] proposed a potential energy method to calculate the effective mesh stiffness and showed that the energy can be divided in three parts: Hertzian, bending and axial compressive energy. These energies can then be used for the calculations of Hertzian contact stiffness, bending stiffness and axial compressive stiffness, respectively. Tian [21] refined this model by adding the shear energy as the fourth part of the potential energy. This analytically obtained time-varying mesh stiffness can represent the effect of changes in the number of teeth pairs in contact and the contact positions between teeth of engaged gears. He also proposed a method to determine the mesh stiffness for different sizes of the crack on the root of a gear tooth. Similar effects on the mesh stiffness can be produced by tooth profile modifications, as highlighted by Chen and Shao [22]. Finally, Iglesias et al. [23] proposed an advanced model for calculation of internal and external gears meshing forces in spur gear planetary. It is then possible to introduce in the model various forms of defects, simply by a proper definition of the mesh stiffness.
The meshing frequency in a gearbox is the frequency at which gear teeth mate together. When the ring gear of a planetary gearbox is fixed, for every complete revolution of the carrier a planet tooth meshes N r times with the ring gear teeth (N r is the ring's number of teeth). Therefore, the meshing frequency can be calculated as follows: where x c denotes the carrier angular velocity (rpm). Furthermore, the angular displacement of the planet gear in one mesh period h m can be calculated by the following equation where T m is the meshing period. In a planetary gearbox several gears with different number of teeth are in contact at the same time (see Fig. 1). All planets are assumed to be identical, thus the behavior and periodicity of each ring-planet or sun-planet meshing are similar although it should be considered these are in different phases with each other and it is also true that each planet's contacts with the ring and sun gears are dissimilar in phase [24]. The contact stiffness related to the generic nth planet is where c rn and c sn are the relative phases between the nth ring-planet and sun-planet meshes, respectively. In addition, c rs represents relative phase between the ring-planet and sun-planet meshes which is identical regardless of which planet is considered [24]. An example of mesh stiffness is given in Sect. 5.1.

Model of the bearings
The outer race, inner race, cage and rolling elements are the key components of a bearing. Figure 2 presents a sketch of the multi-body nonlinear dynamic model that will be used to simulate the vibration response of the planetary gearbox presented in Sect. 2.1. This model was originally developed by Refs. [25,26] which consider two degrees of freedom for the inner race with a fixed outer race. However, in the present model ( Fig. 2) two extra degrees of freedom for outer race are also assumed; as a result, four degrees of freedom of the model comprise the inner raceway displacements x i and y i and the outer raceway  [27] introduce another DOF and they tune its parameters so to represent an high frequency behaviour of the bearings, in accordance with their experimental results. This DOF would add a new equation for each bearing but would not alter the structure of the equations of motion given in Sect. 4. The particular aspect of the flexibility of the bearings, which generates a well separated resonance in the spectra as reported in [27], is not the focus of the present study and has been neglected. The rolling element, the element diameter D b , the pitch diameter D p and the constant operating contact angle n are depicted in Fig. 2a.
In this presented model the rolling elements are supposed massless; thus, the centrifugal forces acting on the balls/rollers are negligible. The flexibility of balls/rollers is modelled by circumferentially distributed radial springs with stiffness k b (Fig. 2b) and the radial clearance is considered as well.
The inner and outer race contact deformations can be combined to calculate the overall contact deformation for the jth rolling element d j as follow where N b is the number of rolling elements, c is the radial clearance and / j is the time variant angular position of the center of the rotating elements. Assuming no slippage or sliding between the components of the bearing, the angular position of each ball, / j , may be calculated based on the races angular rotation and the initial angular position of the first element with respect to x-axis, / 0 .
The cage angular speed can be calculated by using the inner race angular speed, x i , the outer race angular speed, x o , and the geomertry of the bearing [28].
Due to existence of Hertizian contact between balls/ rollers and inner and outer races the reaction force of a roller in position / j is nonlinear and is given by the following expression: Harris [29] suggested exponent n ¼ 1:5 for ball bearings and 1.1 for roller bearings. Springs (balls/rollers) forces will exist only if relative motion between inner and outer race causes moving elements to be compressed. Depending on the sign and amount of clearance it is imaginable that in some cases all balls/rollers will not be in contact simultaneously; therefore, the load zone changes according to the races relative displacements. The contact condition is determined by the positions of the rolling elements. The overall bearing force and stiffness varies with respect to the angular position of balls/rollers and then the total forces in x and y directions can be described as summations of each ball/roller force where c j determines whether the jth ball/roller is in contact, according to

Model of localized faults
When localized faults in inner/outer races or rolling elements are introduced in a bearing model, Eq. (4) should be written in the following form [27] For the inner/outer race defect case, C d in this formulation represents the depth of spalls which is considered to be a function of ball's location inside the defected area. b j determines whether the jth ball/roller is inside defected zone (b j ¼ 1) or not (b j ¼ 0). A rectangular profile is often used for the shape the of spalls but its sharp borders produce large impulsive forces which cause the vibration response of the system to increase too abruptly. To avoid these impulses, a more realistic profile can be used do define the spalls, as pointed out by Liu et al. [30]; in particular, in this paper a Tukey window is generated to model C d (Fig. 2b). The maximum height of this window is chosen equal to the maximum of defect depth on the bearing's inner or outer races. When the ball's diameter is larger than the width of the spall, the maximum depth h theoretically reached by the ball may be calculated as follow (it is an average value for both races) The angular extent of the spall (D/ d ) and position of the spall on the races are also needed in this model to thoroughly define the spall. b j equal to one in Eq. (11) indicates that a rolling element is inside the defected region, thus the compression of the ball and the corresponding force will reduce. As it was mentioned earlier, this variation of force will be revealed by a sudden modification in the acceleration response of the system that can eventually reveal defects. In case of spalls with high depth the ball can even totally lose its contact although this event doesn't happen so often in practice. For the roller/ball defect case, the spall angular speed is identical to the spin speed of a rolling element which can be calculated as follows In this case also a Tukey window is chosen for the spall profile C d (Fig. 2c) and its depth varies as the rolling elements spins. b j ¼ 1 indicates contact between the defect on a roller/ball and the inner or outer race and when b j ¼ 0 there is no contact. For a full rotation of a defected roller/ball, the fault will be in touch with both outer and inner races. Note that in this case the contact duration between the spall and the inner race is longer than with the outer race due to the difference in curvature of the two races [27]. Depending on the race, the total angular contacts between the fault and outer and inner races are where D/ b is the angular width of the defect on the roller/ball (Fig. 2c).
The so far discussed bearing model is employed for sun, planets and carrier bearings. The sun shaft is assumed to be the input, the output is the carrier shaft: the rotational speed of the inner races of the sun and carrier bearings are then equal to their shafts speed. Also, the outer races of sun gear and carrier bearings are considered fixed and have no displacements. Additionally, planets outer raceways have the same displacements as the planets and their inner races are assumed be attached to the carrier. It should also be stressed that the axial movements of the inner and outer races may have an influence on the contact deformation of the rollers and the races. The issue has not been addressed in this paper because of the assumed 2D model (see Fig. 1).

Bearing defect frequency
Bearing defect frequency is the frequency at which rolling elements pass an imperfection on the inner/ outer race or the rolling elements defect approches the raceways. General equations given by Howard [28] are used: where f o is the outer race frequency, f i is the inner race frequency, f bpfi is the inner race defect frequency (BPFI), f bpfo is the outer race defect frequency (BPFO) and f bsf is ball or roller spin frequency (BSF).

Gear defect frequency
Every time a defected tooth of sun gear meshes with a planet tooth a sudden variation will be introduced into the system vibration signal. During one revolution of the sun gear, this faulty tooth will engage with all the planets and therefore the characteristic frequency of sun gear with local fault on a single tooth is calculated as follows where N s is the number of teeth of the sun. When the planet gear has a tooth with a localized fault every time it meshes with ring or sun, a sudden variation will be introduced to the system vibration signal. During one revolution of the planet gear this faulty tooth will engage with ring or sun so this signal modulation occurs once. During one full rotation of the planet the defected tooth meshes twice (with the sun and ring gear) but only one of these contacts may be considered perfect, depending on the damaged tooth side. Therefore, the characteristic frequency of planet gear with local fault on a single tooth is calculated as follows

Equations of motion
The equations of motion may be written based on the mentioned lumped parameter model of a planetary gearbox and the bearings model. It is worth noting that gyroscopic and centrifugal forces are also considered in this dynamic model. The sun equations of motion are: where T i is the input torque of the system and F spn represents the gear mesh force between the n-th planet and sun gears where a is the pressure angle of gears and N is the number of planets. The planets are assumed equally spaced and W n represents the angular distance between the planets. F sbx and F sby represent the sun bearing force in x s and y s directions based on Eqs. (8) and (9) and / sj , the angular position of sun bearing balls, is calculated with Eq. (5).
The ring equations of motion are: m r € y r þ c ry _ y r þ k ry y r þ X N n¼1 F rpn sinW rn where F rpn represents the gear mesh force between the n-th planet and ring gears The planets equations of motion (n ¼ 1; . . .; N) are: where F cpxn and F cpyn represent the planets bearing force in x pn and y pn directions based on Eqs. (8) and (9) and / pj , the angular position of the planet-bearing balls, is calculated with Eq. (5).
The carrier equations of motion are: where T o is the output torque of the system. F cbx and F cby represent the carrier bearing force in x s and y s directions based on Eqs. (8) and (9) and / cj , the angular position of the carrier bearing balls, is calculated with Eq. (5).
It is stressed that the equations of motion are written in the rotating frame of reference Oxy and their results are therefore almost impossible to measure in an experimental rig. In practice, vibration signals of planetary gearboxes can in fact be collected via accelerometers mounted on the casings, being in most cases the internal parts not accessible. Once the previous equations have been solved, it is necessary to project the results in the fixed frame of reference OXY (see Fig. 1). In particular, the results discussed in the following section arise from the accelerations of a fixed point on the ring gear, to simulate an actual condition. Gravity can be added to the equations of motion but it has been numerically verified that the presented results would not be modified. In fact, with the parameters set in Table 1, the effect of the weight force is negligible in comparison to gear meshing forces and, as a result, does not influence the computed time histories, i.e. the accelerations of the output point.
To the best knowledge of the authors no published work, except [18], has been found in the literature which try to simulate the effect of bearing localised faults in a planetary gearbox. The present model has indeed some differences with may be worth noting. First, Jain [18], takes into account linear and timeinvariant mesh stiffness between the gears, which allows him to define time-invariant mass and stiffness matrices. Consequently, he is able to define and solve an eigenvalue problem and to eventually determine the frequency response functions (FRFs) of the system in terms of modal properties. On the contrary, our proposed model uses time-variant forces and also considers mesh phasing, which can have a substantial effect on the result [24]. A constant mesh stiffness also makes it impossible to introduce any defect on gears, since they are modelled by variation on mesh stiffness (see Fig. 3). Moreover the nonlinearity of bearings contacts in addition to fluctuation of the bearing total stiffness, as a result of the cage rotation, are considered. The achieved equations are then non-linear and time-variant and can't but numerically be solved in the time domain. It also may be observed that computing a FRF in terms of receptance is not too realistic because it is impossible to isolate a single force, as required by the definition of the FRFs. Second, the signature of bearing localised faults is analytically modelled in [18] as a modification of the bearing force, which limits the ability to model bearing faults with different shapes. This approach is here replaced by the model in the Sect. 2.4, which describes the damage in terms of variation of the contact pattern. Third, the size of the bearing load zone is considered constant in [18] which is not the case here, since it differs for bearings with various number of rolling elements. Finally, the flexibility of the ring gear is modelled in [18] by using a modal expansion, which has here been ignored. The plan is to expand the model so to connect the ring gear, the bearing of the sun gear and the bearing of the carrier to a flexible structure but, at this stage, the elasticity of the ring gear has been disregarded.

Results and discussion
In this section, the set of nonlinear and time-variant equations derived in Sect. 4 is solved to obtain vibration signals of each gear. In this regards, a code in the Wolfram Mathematica software environment has been developed to numerically solve this system of ordinary differential equations (ODE) by means of the built-in NDSolve function. The gears and bearings parameters are included in Tables 1 and 2 respectively.

Gears with cracked teeth
The approach of Ref. [31] is used to analytically evaluate the time-varying mesh stiffness of externalexternal and external-internal gears of our planetary gear set. Figure 3 shows the mesh stiffness of perfect and cracked (10 percent of tooth root is cracked) teeth in the sun-planet and planet-ring coupling for a meshing duration corresponding to the average number of gear teeth pairs in contact while a tooth comes and goes out of contact (contact ratio). As Fig. 3 demonstrates, the amplitude of mesh stiffness waveform decreases for damaged gear teeth due to the thickness reduction of the cracked tooth. In this research c rs is 0.5 and c sn and c rn for n ¼ 1; 2; 3 are 0, 2/3, 1/3 and 0, -1/3, -2/3, respectively (see Sect. 2.2). The mathematical description of gears contact damping coefficient is complex. In this study a simplified damping model, Ref. [32], is assumed to determine the effective damping factor as follows: where f is the contact damping ratio. The value of f is between 0.03 to 0.17 [32] and in this study an average value of 0.10 is selected. Figure 4 illustrates the acceleration signal of the ring gear (in the fixed frame OXY) when the crack is seeded in a single tooth of the sun gear (Fig. 4a) or planet gear (Fig. 4b). When tooth cracks are present, impulsive signals can be observed in time domain. The time duration between two consecutive disturbances in the acceleration signal is equal to the gear defect period. Based on the Eqs. (18) and (19) sun and planet gear defect frequencies can be calculated as 18.46 Hz (0.054 s) and 8.42 Hz (0.119 s) respectively.

Defected planet-bearing
In this section effects of faulty inner race, outer race and rolling element of the planet-bearing will be mainly considered.
In a planetary gearbox forces acting on carrier through planets bearings or on sun gear through the planet-sun gear mesh counterbalance each other and the resultant forces on carrier/sun bearing is negligible. On the other hand, radial load on planets bearings, which transmit the torque are much higher. As a result, in contrast to sun and carrier, planets bearings exhibit a high failure rate and are considered as one of the most critical components in planetary gearboxes. Therefore, the focus of the paper is on defects of planets bearings in the gearbox.

Fault on inner race
For the healthy gearbox, the acceleration of the ring gear in X direction for zero clearance is shown in Fig. 5a. Also, Fig. 5b displays the ring acceleration result when a 4-degree defect (D/ d ¼ 4 o ) is seeded in the inner race of the first planet bearing. The localized fault is positioned along the rotating y axis of the planet-bearing and does not move with respect to carrier and bearing load-zone.
For the planets bearings the inner race frequency is equal to carrier frequency (f i ¼f c Þ and the outer race is fixed to the planet (f o ¼f p ¼ À2:37f c Þ. Therefore, the ball pass frequency of planet inner race (BPFI) when the carrier rotational speed is 150 rpm can be calculated according to Eq. (15) as 41.39 Hz. Dashed Contact angle, n 0 Fig. 5b represents the time duration in which a roller goes through the spall and the spacing between every two successive occurrences is equal to defect period (T i ¼1=f bpfi ¼ 0:024 s) which also verifies the validity of the implemented model. As the fault is always inside the load zone its effect is present in every passage of balls over the fault. Although the bearing fault signature can be presumed inside the dashed box, the signal is mainly dominated by gears meshing components. To attain clearer perception of bearing defect signature, zoomed portions of Fig. 5a, b are plotted in Fig. 5c. Disturbance of the acceleration signal which is due to the defect on the planet-bearing can be seen more evidently.

Fault on outer race
Planets are responsible to transmit the input torque from sun gear to carrier and therefore forces on their bearings in y direction, which causes the carrier to rotate, are much larger than forces in x direction. When positive or zero clearance exists in planets bearings there will be no sign of defect once the fault is out of the load zone. In this simulation, to represent the effect of outer race spall on output signal more clearly, the clearance is set to an arbitrary value of À5 lm in the planet-bearing model since the negative value of clearance or preload could be generated by elastohydrodynamic lubrication films (EHL) [27].
Ring gear acceleration signal in X direction in case of healthy and faulty planet-bearing, 4-degree rotating spall, are displayed in Fig. 6a, b respectively. To clarify the influence of outer race spall in the ring gear signal, zoomed portion of Fig. 6a, b for the period in which one of the balls is inside the defected area are simultaneously represented in Fig. 6c. By using Eq. (16) defect frequency for the planet outer race (BPFO) is calculated as 25.97 Hz. Defect time period in this case is equal to T o ¼1=f bpfo ¼ 0:039 s which is shown in Fig. 6b as time difference between two successive defected areas of acceleration signal For the case of planet-bearing with a defected ball, BSF is calculated by using Eq. (17) as 34.87 Hz. Since the balls defects have somehow the same characteristics as the outer race defects, the time domain results are not depicted to not be repetitive.

Frequency analysis of defect signals
This subsection is devoted to the analysis of the frequency spectrum of signals associated to planetbearing defects. The objective is to highlight the signal components generated by damage and their sources. The undamaged bearing is considered to generate a reference signal. Hence, the difference between the To determine the frequency spectrum, a Fourier transform is performed on the residual signals. Figure 7 shows the spectra of residual vibration response of the ring gear due to the presence of a single defect on the inner race of the planet-bearing.
Inner race fault frequency, f bpfi , and its harmonics can be seen in Fig. 7a. Each of them encompasses a cluster of sidebands separated by the carrier rotation frequency, f c : a magnified part of Fig. 7a is represented in Fig. 7b. These sidebands are produced as a result of carrier rotation. In this case, as mentioned earlier, the ring gear is fixed so the centers of planets revolve along with the carrier. This revolution causes the transmission path between the planet-bearing and the signal acquisition point on the ring gear to vary with time. Due to this variation, the amplitude of planet's bearing vibration signal is modulated and generates the sidebands around the damage frequency and its harmonics. It must be addressed that a spectrum with proper frequency resolution is needed to detect the damage frequencies and their sidebands, so time duration of calculated (or measured) data should be chosen long enough to generate a reliable frequency spectrum.
Frequency spectrum for planet-bearing with defected outer race and ball are shown in Figs. 8 and 9 respectively. Similar to the previous case, a cluster of sidebands is present around the outer race damage frequency, f bpfo , and ball damage frequency, f bsf , and their harmonics. But as it is shown in Figs. 8b and 9b more sidebands around damage frequencies exist in comparison to the inner race fault since in these cases two frequencies are responsible for creation of these sidebands. As in the former circumstance, the carrier rotation changes the transmission path and is the first source of planet-bearing signal amplitude modulation so that f bpfo and f bsf are modulated by f c . When a spall is located on planet-bearing inner race, its position doesn't change relative to rotating frame xy because both are attached to the carrier. However outer race fault is fixed to the planet gear and relative rotation occurs between it and the carrier; moreover, relative rotation between ball and carrier occurs as a result of the cage rotation. Consequently, the following events occur.
First, unlike the inner race fault, the sudden force which is caused by the defective planet-bearing outer race or ball will rotate relative to the xy reference frame. The amount of these relative rotational speeds is equal to x cage in case of ball defect and can be calculated by subtracting the speed of the planet (x p ) from the speed of the carrier (x c ) in case of outer race defect.
Second, if the planet-bearing clearance is zero/positive or negative but not sufficient to prevent the ball/ roller from losing its contact with the races while it passes through the fault or when the ball defect is in contact with the races, the magnitude of force generated by the defect is different according to the position of the spall relative to the load zone (i.e. amplitude of force will be non-zero inside the load zone and zero when the spall is out of the load zone).
Third, planets vibration signals are transmitted to the ring gear via ring-planet engaging teeth. Due to the spall rotation along with the planet gear or bearing cage the angle between the defect force and ringplanet mesh varies. Maximum force will be transferred to the ring gear when the defect force and the gear meshing spring are in the same direction. On the other hand, perpendicularity minimizes the transmitted force [18].
According to above mentioned phenomena, the rotation of defect together with the outer race or cage changes magnitude and direction of its force. Consequently, these variations generate a second frequency of modulation which is equal to f pc ¼ f p À f c for the outer race spall and f cage for the ball spall. As two modulation frequencies exist for the outer race and ball fault cases, the cluster of peaks around each damage frequency (sidebands), which are shown in Figs race fault or f c and f cage for the ball fault and can be calculated as follow [18] where It is worth noting that as bearings balls move, their total stiffness varies. In case of planet-bearing, frequency of this variation is equal to inner race ball passage frequency. Thus if the preload and number of balls are not properly selected, even for flawless bearings, vibrations occur at this frequency but increasing the number of balls decreases the amplitudes of bearing vibrations and ball pass frequency (BPF) will reduce accordingly.
The damage frequencies and sidebands presented for both planet inner race, outer race and balls faults in this section comply with the theoretical and experimental finding of Ref. [18]. This agreement also verifies the reliability of the bearings and gearbox models as well as accuracy of the simulation results. Figure 10 shows the frequency spectrum of the ring gear acceleration signal for healthy and defected planet inner race bearings. In Fig. 10b, c the sizes of defects are 2 and 4 degrees respectively and Fig. 10a represent the frequency spectrum of the healthy gearbox. Around the meshing frequency (160 Hz) sidebands emerge when a defect is present. Although the amplitudes of meshing frequency and its harmonics are barely affected by the bearings faults, amplitudes of sidebands increase as the size of spall grows. These collection of frequencies around the meshing frequencies are combination of BPFI and its sidebands as discussed in the previous subsection. The spectrum of the gearbox signal with healthy, 2 and 4 degrees faults on the outer race of one planetbearing are also illustrated in Fig. 11. Notice that Fig. 10a differs from Fig. 11a because of the presence of clearance. In this case, the spectra comprise more sidebands as explained in Sect. 5.2.3. The BPFO frequency and its sidebands appears in vicinity of the meshing frequencies.

Frequency spectrum analysis of the ring acceleration signals
The spectrum for healthy and defected planetbearing balls are depicted in Fig. 12. The maximum balls spalls depths are chosen equal to maximum depths of races defects, h, to have comparable results. The BSF frequency and its sideband are also present in the spectrum. These spectra are similar to those generated by outer race defects case because in both circumstances defects are in rotation relative to the carrier and therefore two sources of modulation exist as explained in the previous subsections.
In all above cases, the amplitudes of meshing frequency and its harmonic are higher than the amplitude of the bearing defect frequency and its harmonics. This is caused by the fact that planetary gearboxes transfer a large amount of torque from sun gear to carrier through their gear meshes causing rather high energy gear mesh frequencies in comparison to the energy level of signals generated by bearing defects.

Condition monitoring of the gearbox
Diagnosis of defects on planets bearings at their early stages has always been problematic. Many statistical features have been developed for condition monitoring of gearboxes and have been frequently used as an indicator of bearings and gears faults presence and growth. Therefore, the vibration signal from the simulation is processed to calculate a set of statistical features. The accelerations for different levels of planet-bearing defects on inner races, outer races and rolling elements will be evaluated by implementing the root mean square (RMS), kurtosis and M8A features to investigate the effectiveness of statistical indicators in detection of planets bearings defects. Root mean square (RMS) might be one of the most commonly used indicators in vibration monitoring. It is defined as where x i is the data sample, N and x are the length and average of x respectively. The kurtosis is the fourth normalized moment of a signal and provides a measure of its peakedness. It is given by The parameter M8A uses the eighth moment normalized by the variance to the fourth power and is given as Many research work has been carried out in the field of statistical fault detection and diagnosis. Typically time domain signals are the inputs of statistical indicators [33]. In this work these methods are applied on the ring acceleration in time domain but when the sizes of defects are small these indicators are not reliable and cannot detect the presence of the fault because the signal is highly dominated by the gear components.
As it was discussed earlier, frequency sidebands appear about gear mesh frequencies when a bearing fault is introduced and amplitudes of sidebands intensify as the spall grows. These variations are believed to be worthwhile for fault diagnosis of planetary gearboxes. Therefore, the condition indicators on the ring acceleration are exploited in frequency domain to investigate the planet-bearing defect. The procedure is the same as the time domain analysis except that data from Fourier transform of the ring gear time signals are take into account. For every result discussed in this section, white noise (with RMS  Fig. 12 Frequency spectrum a healthy gearbox b defected planet balls, equal to 15 percent of the acceleration RMS) has been added to the original time sequence to test the performance of these indicators when background noise is not negligible.
The calculated condition indicators for 6 cases, no defect and 5 different sizes of the spalls (D/ d ¼ 1, 2, 4, 8 and 15 degrees), are calculated. For the defected balls D/ b ¼ 45 o and the maximum depth of the ball defect is selected equivalent to maximum depth of the races spall h. For defects on the outer race, inner race and balls of planet-bearing, the difference between the indicators in defected and healthy cases are calculated then these differences are divided by the healthy bearing indicator. The result are plotted in Fig. 13.
As it is shown in the Figs. 10, 11 and 12, when the spall advances not only new peaks appear in the spectrum but also the amplitudes around meshing frequencies ascend generally. Therefore, all the statistical indicators quantities rise with the growth of the planet's defect size (Fig. 13). As the average of the data increases the value of RMS gets larger correspondingly, while kurtosis and M8A provide measures of the peakedness of data sequences. A higher positive value indicates a peaked distribution while a lower negative value indicates a flat distribution. When the fault is present, as more peaks appear around the gear mesh frequencies the values of kurtosis and M8A decrease i.e. the effect of gear mesh frequencies deviation on variance reduces because of sidebands peaks. Therefore, the absolute values of the indicators are presented in Fig. 13. The M8A indicator is more sensitive to the presence of the defects and performs better on faults detection, with a potential to diagnose spalls even in their emerging stages. Another observation is that the lower frequencies are more sensitive to the presence of the faults. Lower natural frequencies of the system have larger amplitudes therefore they more effectively amplify the amplitudes of the meshing and defect frequencies. Moreover, as the frequencies amplitudes escalate around the natural frequencies of the system the effect of noise becomes less substantial. At the end of this section it is worth noting that increased sensitivity may not always be a desired characteristic because an excessively sensitive parameter will yield many false alarms. For this reason, the most useful damage indicator is not definitely the most sensitive one.

Conclusion
In this paper interaction between gears and bearings of a planetary gear set in the presence of faults is investigated by developing a lumped parameter model. The capability of simultaneously simulating the response of the system for different sizes, locations and profiles of bearing and gears defects is the significant advantage of this proposed model. The acceleration signal of ring gear for healthy and defected bearings and gears are calculated and compared. Moreover, frequency components of signals associated to planets bearings defects for inner/outer race and rolling elements are discussed. In case of inner race faults, the orbital rotation of planets and in case of outer race/roller faults rotation of spalls relative to carrier as well as planet orbital rotation are sources of amplitude modulation of gearbox components signals and as a result sideband clusters around the bearing damage frequencies. Furthermore, effects of faults on planets bearings and growth on frequency spectrum of the ring gear acceleration signal are investigated. Influence of different levels of inner/outer race and rollers planet-bearing defects on dynamic response of the system are evaluated by implementing statistical indicators (RMS, kurtosis, M8A). It is concluded that these indicators are more effective when they are applied to the frequency domain data rather than to time domain signals. The quantities of these three statistical indicators increase as the spalls grow. M8A is the most sensitive feature and seems to have the potential to detects damages in the initial stages.
For future work, results of the reported model in this article can be used to test damage detection techniques especially when different gears and bearings defects are present simultaneously. Furthermore, this dynamic simulation can be modified to take into account flexibility of gearboxes casings, slippage in bearings and also variable speeds of gears.