Symmetry breaking and chaos-induced imbalance in planetary gears

The goal of the present paper was a complete analysis of the dynamic scenario of planetary gears. A lumped mass two-dimensional model is adopted; the model takes into account: time-varying stiffness; nonsmooth nonlinearity due to the backlash, i.e., teeth contact loosing; and bearing compliance. The time-varying meshing stiffness is evaluated by means of a nonlinear finite element model, which allows an accurate evaluation of global and local teeth deformation. The dynamic model is validated by comparisons with the most authoritative literature: linear natural frequencies and nonlinear response. The dynamic scenario is analyzed over a reasonable engineering range in terms of rotation speed and torque. The classical amplitude–frequency diagrams are accompanied by bifurcation diagrams, and for specific regimes, the spectral and topological properties of the response are discussed. Periodic, quasiperiodic and chaotic regimes are found; nonsmooth bifurcations lead period one to period two trajectories. It is found that the bearing compliance can influence the natural frequencies combination magnifying the modal interactions due to internal resonances and greatly enlarging the chaotic regions. It is evidenced that the chaotic response indices a symmetry breaking in the dynamical systems. The physical consequence is that the planetary gearbox under investigation, which is perfectly balanced for each position, can suffer of a big dynamic imbalance when chaotic regimes take place; such imbalance gives rise to alternate and unexpected high-level stresses on bearings.

are found; nonsmooth bifurcations lead period one to period two trajectories. It is found that the bearing compliance can influence the natural frequencies combination magnifying the modal interactions due to internal resonances and greatly enlarging the chaotic regions. It is evidenced that the chaotic response indices a symmetry breaking in the dynamical systems. The physical consequence is that the planetary gearbox under investigation, which is perfectly balanced for each position, can suffer of a big dynamic imbalance when chaotic regimes take place; such imbalance gives rise to alternate and unexpected high-level stresses on bearings.

List of symbols
x, y x , y-translational coordinate (µm) Θ Rotational coordinate (rad) k s Translational sun bearing stiffness (equal for x, y) (N/m) Ψ n Angular position of the nth-planet with respect to the x-axis (rad) α s Working pressure angle for sun-planet mesh (rad) α r Working pressure angle for planet-ring mesh (rad) r bs , r bn , r br Base radii of (sun, nth-planet, ring) (m) k sn Mesh stiffness of sun-nth-planet (N/m) k su Torsional sun bearing stiffness (N m/rad) k r Translational ring bearing stiffness (equal for x, y) (N/m) k rn Mesh stiffness of nth-planet-ring (N/m) r c Carrier radius (sun-planet center distance) (m) k p Translational planet bearing stiffness (equal for x, y) (N/m) k pu Rotational planet bearing stiffness (in this paper k pu = 0) (N m/rad) k c Translational carrier bearing stiffness (equal for x, y) (N/m) k cu Rotational

Introduction
Planetary or epicyclic gear trains are widely used in many automotive, aerospace and marine applications; they are effective power transmission systems when high torque-to-weight ratios, large speed reductions in compact volumes, coaxial shaft arrangements, high reliability and superior efficiency are required [1]. Gear vibrations are primary concerns in most planetary gear transmission applications, where the manifest problem may be noise or dynamic forces. The most important source of vibration in planetary gears is the parametric excitation due to the periodically timevarying mesh stiffness of each sun-planet and ringplanet gear, because the number of tooth pairs in contact changes during gear rotation. This mesh stiffness variation parametrically excites the planetary gear system, causing severe vibrations when a harmonic component approaches one of the natural frequencies (or their linear combinations). Under certain near-resonant operating conditions, gear systems can experience a teeth separation that induces nonlinear effects such as jump phenomena and subharmonic and superharmonic resonances with dramatic effects on the dynamic response [2]. These phenomena have been deeply investigated in geared systems during the last 20 years [3][4][5][6][7]. Cunliffe et al. [8], Botman [9] and Kahraman [10][11][12] presented models to estimate the natural frequencies, vibration modes and dynamic forces of planetary gears. Ambarisha and Parker [1] validated the effectiveness of a lumped-parameter model to simulate the dynamics of planetary gears. In Ref. [1], responses from the dynamic analysis using analytical and finite element models are successfully compared qualitatively and quantitatively.
Parker and Lin [19] found that the mesh phasing has a dramatic influence on the static and dynamic behavior of planetary and epicyclic gears. Linear system analysis helps to understand resonances and parametric instabilities due to periodically varying gear mesh stiffness [20], mesh phasing [19,20] and properties of vibration modes [14,15,23,24].
Some researchers considered two-stage planetary systems [24,25]. Guo and Parker [26] derived important relations among the relative phases between any two gear meshes in a compound planetary gear.
Al-shyyab and Kahraman [27] solved the nonlinear equations of planetary gear motion semi-analytically using multi-term harmonic balance method (HBM) in conjunction with inverse discrete Fourier transform and Newton-Rapson method. The nonlinear dynamics of planetary gears involving teeth wedging and bearing clearance was investigated by Guo and Parker [28]. Sun and Hu [29] investigated the nonlinear dynamics of a planetary gear system with multiple clearances; the theoretical results from HBM are verified by using the numerical integration. When the planetary gearbox operates under some inappropriate conditions, such as inadequate lubrication, poor specifications, material defects and improper manufacturing or installation, it is likely to cause gear teeth faults. Once the tooth faults appear, the dynamic performances of the gearbox are bound to be affected, with undesirable changes in the dynamic behavior and serious reduction in the service life of the planetary system [30][31][32][33]. Recently, Gu and Velex developed a dynamic model including three-dimensional effects and errors [34,35], and Li et al. [36] proposed a deep investigation on chaotic effects in planetary gears.
This paper presents a dynamic model to simulate the dynamic behavior of a single-stage planetary gear system with time-varying mesh stiffness and backlash. The complex dynamic scenario of a three-planet gearbox is investigated in detail. A bifurcation analysis is performed to explore the dynamic scenario (periodic, quasiperiodic and chaotic), with a special attention to symmetry breaking phenomena that are extremely interesting in planetary gears as they can cause additional imbalance-induced stresses. Numerical analyses are carried out over meaningful mesh frequency ranges. The analysis is completed with time histories, spectra, phase portraits and Poincaré maps of the most interesting regimes.

Dynamical model
The physical model of a single-stage planetary gear set is shown in Fig. 1. The system is made of four types of elements: sun gear, ring gear, N planets and carrier. Here, the modeling is plane, i.e., each element has three degrees of freedom: two displacements and a rotation. The centers of the different elements of the system are free to move in the plane, and each component has translational and rotational degrees of freedom. The total number of degrees of freedom is (3N + 9); the model includes time variation of gear mesh stiffness (depending on the reciprocal angular position of two meshing gears) k sn and k rn (n = 1, 2, . . . , N ), back-lash nonlinearities of mating gears and bearing compliance (no clearance is considered for bearings).

A single-stage planetary gear set
In Fig. 1, the planar physical model for the planetary gearbox vibration analysis is presented. Each planet is mounted on a rigid carrier through a flexible bearing; sun, ring and carrier bearings are connected to the ground, while the planet bearings are connected to the carrier. The planets are free to rotate with respect to the carrier; the sun is the system power input, and it is also free to rotate; the rotational degree of freedom of the ring is locked; the rotation of the carrier is the system power output, and it is constrained by a rotational bearing stiffness. The gear-shaft bodies, the ring and the carrier are assumed rigid; the compliant elements of the meshing gear teeth and bearings are represented by suitable nonlinear time-varying springs. The origin of the global coordinate system is at the undeflected position of the center of the sun. Positive planet position angles ψ i (i = 1, 2, . . . , N ) are measured counterclockwise from the first planet, i.e., ψ 1 = 0, see Fig. 1.

Details of the single-stage planetary gear set model
The gear elements, S, R and C, are constrained to the ground by means of translational and torsional linear springs of stiffness magnitude (k s , k su ), (k r , k ru ) and (k c , k cu ), respectively. Each spring is associated with a linear viscous damper, proportional to the mean value of stiffness and inertia, in order to simulate energy dissipation of the system. The planets are connected to the carrier by means of translational stiffness and damping. A proper choice of the magnitudes of the stiffness constraints allows to represent any different power flow scheme. Each gear body i(i = S, R, C, P 1 , . . . , P N ) is modeled as a rigid disc of mass M i , mass moment of inertia I i , base radius r bi and torsional displacement θ i . Planets are located at radius r c angular positions ψ i ; here, the planets are considered equally spaced circumferentially. Excitation torque T s is applied to the sun gear. The mesh of sun (or ring) gear with a planet P i is represented by a periodically time-varying stiffness element k sn (t) (or k rn (t)) subjected to a piecewise linear backlash function f s (or f r ) that includes a clearance of amplitude 2b s (or 2b r ). The vibrational These coefficients are proportional to the mean stiffness coefficients and also to the mass of gears such that C = αM + βK, where C, M and K are damping, mass and stiffness matrices, respectively, and α and β are coefficients given by natural frequencies and modal damping ratios of the system, such that the damping ratio is ρ i = 1/2(α/ω i +βω i ) for the generic ith mode.

Equations of motion
The basic dynamical equilibrium equations contain (3N + 9) nonlinear ordinary differential equations, where N is the number of planets; e.g., when N = 3 they will be 18 coupled equations. The equations of motion of the model shown in Fig. 1 are written using Newton-Euler equations, and they are placed in canonical form.

Sun gear equations
+ k sn f sy × cos (ψ n −α s ) × r bs In these equations, x and y and θ are the translational and rotational degrees of freedom, respectively. T s is the constant input torque applied to the sun gear and f sx , f sy and f sθ are the piecewise linear displacement functions for sun-planets meshing, which are defined as follows: where, When | sn | < b s , the teeth are separated, and when sn ≤ −b s , backside contact occurs (Fig. 2).

Ring gear equations
+ k rn f r y × sin (ψ n + α r ) × cos(ψ n + α r ) where, f r x , f r y and f r θ are the following: Similar to the sun-planet contact, tooth separation occurs when | rn | < b r , while if rn ≤ −b r the contact is on the opposite flank (Fig. 2).

Planets gear equations
(8c) In Eqs. (1)(2)(3)(4)(5)(6)(7)(8), k sn (t) and k rn (t) are the mesh stiffnesses between sun-planet and ring-planet gears, respectively; they are periodically time varying at the mesh frequency. As the number of teeth pairs in contact within the mesh cycle changes with the rotation of the gear, the mesh stiffness varies accordingly. The piecewise linear spring generates nonzero mesh force only for positive relative mesh deflection (teeth in contact). The combination of the static torque and time-varying stiffness gives rise to the excitation of the system, and the periodic stiffness has the fundamental frequency ω m = ω s z s z r /(z s + z r ), where ω s is the angular velocity of the sun gear and z s and z r are the teeth numbers of the sun and the ring gears, respectively. Periodically time-varying mesh stiffnesses are written in Fourier series form: where h s and h r are the number of harmonic terms used to describe the periodic functions k sn (t) and k rn (t); k r p are the harmonic coefficients of Fourier series. k (1) sp and k (1) r p are the mean values of the coefficients; γ sn represents the phasing between sun-planets, which means relative phase between nth sun-planet mesh and first sun-planet mesh, i.e., γ s1 = 0; the same definition applies toγ rn , which is relative phase between nth ringplanet mesh and first ring-planet mesh,γ rn = γ sn +γ sr represents the ring-planets phasing; γ sr is the sunring phasing, i.e., the phase angle between the nth ring-planet mesh and the nth sun-planet mesh. For the present study, γ sr is independent from planet number; so, γ sr is phase angle between stiffnesses of meshes sun-first planet and ring-first planet. See Ref. [21] for a comprehensive description on mesh phase relations in planetary gears.

Time-varying mesh stiffness evaluation
The gear meshes are considered as linear springs with backlash (nonlinearity) and periodically varying stiffness, due to the variation of tooth contact conditions and position of contact points on each tooth. This time-varying mesh stiffness is evaluated by Helical-Pair software [11,38]. HelicalPair is a software developed in the Centre INTERMECH MO.RE. (Aster, High Technology Network of the Emilia Romagna Region). Linearized gear mesh stiffness (without considering backlash effects) is evaluated separately both for external gears (sun-planets mesh) and internal gears (ringplanets mesh). HelicalPair acts as a preprocessor and postprocessor for a FE solver, namely MSC Marc software. Figure 3a, b illustrates finite element model produced by HelicalPair software for external and internal gear meshes. Figure 3c, d presents details of gear tooth meshes; backlashes are observable in these figures.
The mesh stiffnesses are dependent on the forcing condition (torque applied). For external gear meshes, the torque T s = T /N is applied to the sun, where T is the sun torque of the planetary system, the transmission error is calculated using HelicalPair (u s = r s θ s ), and the mesh stiffnessk s1 at a particular mesh position is evaluated ask s1 = T s /(r s u s ). This calculation is repeated at multiple steps within a mesh cycle. A similar process is used for the ring-planet mesh; the torque T p = T s z p /z s is applied to the planet. In each case, the center of the driven gear is fixed when the stiffness is evaluated.
Note that the present calculations of the applied torques are valid when the load is uniformly distributed among different planets (balanced gearbox). In actual gearboxes, it is impossible to guarantee uniform loading when more than three planets are present, due to the quality accuracy of workmanship [33].

First validation: pure rotational model
In this section, the results are validated by comparisons with Ref. [2], where a pure rotational nonlinear model was considered. The gear tooth flanks are pure involutes without tooth modification. The parameters of the gearbox are given in Table 1. The gearbox is made of three equally spaced planets, i.e., N = 3 and ψ i = 0, 2π/3, 4π/3. The half backlash values for both meshes are b s = b r = 0.3 mm [2]. A constant sun gear torque of 1,130 Nm is applied. The ring gear is fixed in the present study. The carrier rotational dof θ c is constrained to be zero under the assumption that there is a large output inertia (as usual for many applications), so the rigid body motion is removed [2].
The average stiffness values of the sun-planet and ring-planet meshes are: k   Table 2. The values of the gear mesh damping constants C sn (t) (or C rn (t)) are estimated such that they correspond to the damping ratios given in Table 3.
All sun-planet meshes are in phase with each other, as well as all ring-planet meshes, i.e., γ sn = γ rn = 0. This design can be realized when the summation of sun and ring teeth numbers is an integer multiple of the number of planets [2]. In the absence of chaotic responses, caused by nonlinear effects, this arrangement presents symmetric force distribution, i.e., perfect balancing on the sun (zero forces acting on the bearings).
The natural frequencies of the present model (reduced to pure rotational motion by elimination of translational degrees of freedom) are presented in Table 4. Good agreement between Ref. [2] and the present model is found. The nonlinear model is validated considering regimes where the amplitudes of vibration cause teeth separation, i.e., nonlinear behavior. Figure 4 shows the RMS (root mean square) of sun rotational vibration (multiplied by sun base radius) versus the meshing frequency; such results are obtained by numerical integrating Eqs. (1)(2)(3)(4)(5)(6)(7)(8); the mean value of time histories is removed before calculating the RMS. The algorithm for numerical integration is RADAU5; the number of mesh periods considered for transient response is 950 (initial response removed from the analysis), the number of substeps is 100 per each period, and the number of mesh periods considered for evaluating the RMS is 50. The mesh frequency is dependent on the sun rotational speed by ω m = z s z r /(z s + z r )ω s .
Results for both increasing and decreasing speed sweeps are plotted in Fig. 4, and only the sun rotational dof is represented for the sake of brevity. There are jump phenomena at the first (ω 1 ) and the fourth (ω 4 ) resonances. The resonance peaks lean to the left, implying softening nonlinearity induced by tooth separation. There are additional resonance peaks around 8,000 Hz, 1,800 Hz (the first distinct mode) and below 1,600 Hz. These peaks are the combined effect of parametric instability from higher harmonics of mesh stiffness fluctuation and nonlinear subharmonic and superharmonic resonances of the first and second distinct modes [2].

Second validation: rotational-translational model
The full model for nonlinear dynamics of planetary gears is investigated here using all degrees of freedom for sun, planets, carrier and ring gears. Comparison between the rotational-translational bulk model with the results of Ref. [1] is presented. The natural frequencies of both models are listed in Table 6 (gearbox data are available in Tables 2, 3   The eigenvectors correlated with the natural frequencies of Table 5 are plotted in Fig. 6. These eigenvectors are normalized by unity. Table 7 presents the description of each degree of freedom. The modes Fig. 6 Six first mode shapes of the system corresponding to Table 6. Definition of degrees of freedoms is given in Table 7

Nonlinear vibrations: bifurcation analysis
The nonlinear response of the rotational-translational bulk model of planetary gear set is now analyzed in detail by means of bifurcation diagrams, Poincaré maps, phase portraits, time responses and spectra. The complex dynamics is analyzed via direct simulation. In order to obtain the Poincaré maps, time histories are sampled with the same frequency of the excitation. Poincaré sections are useful tools for analyzing chaotic responses. Moreover, bifurcation diagrams of the Poincaré maps can be drawn by varying one of the parameters, which govern the response of the system (the meshing frequency here). The equations of motion are integrated for a sufficiently long time in order to eliminate the transitory motion, and then a certain number of periods are considered. This procedure is repeated for different values of the parameter, using the final state of the system in the previous analysis as initial condition. The procedure allows a particular, stable solution to be followed when a parameter is varied regularly. Numerical parameters of Tables 2, 3 and 5 are considered. Two different translational stiffnesses of planet bearings are investigated: K p = 2.19 × 10 9 N/m (from Ref. [2]) and K p = 2.19 × 10 8 N/m. Natural frequencies of each case are listed in Table 8.
The comparison between two cases in Table 8 shows that stiffness of the planet bearings has bigger influence on low natural frequencies. The higher natural frequencies are similar in both cases, while the fundamental frequency of the first case is more than two times the one of the second case. The symmetry of the system gives rise to double modes.
The main excitation frequency of the system is the meshing frequency of the planetary gear. The amplitude of external excitation is related to the external torque on the sun gear. Behind the main excitation, there exist internal resonances due to the integer ratios between the natural frequencies of the system. We show the change in nonlinear phenomena between the systems with and without internal resonances.
For the two cases, the internal resonances are summarized below: Case 1 of Table 8: ω 7 ∼ = 2 × ω 4 , ω 15 ∼ = 6 × ω 6 Case 2 of Table 8: In the case 2, the internal resonances involve the fundamental frequency, so that the resulting vibration is more significant.
The bifurcation diagrams are extracted by varying the mesh frequency; 750 frequencies (upward) are analyzed, for each one the integration is carried out for 1,000 periods, the last 50 periods are recorded and the initial 950 periods are disregarded in order to remove the transient vibrations; when the frequency is changed, the new simulation is carried out using the final state of the previous simulation as initial condition with a small perturbation; the sampling frequency is equal to the excitation frequency in order to build the Poincaré maps. Figures 7 and 8 show bifurcation diagrams versus excitation mesh frequency for the cases 1 and 2, respectively. Figures 7a and 8a Fig. 7 (for case 1), it can be observed that the bifurcation diagram crosses a 2T region at ω m ∈ (10,200−10,700) Hz, where the period of response is twice the excitation period; this 2T region disappears for case 2. The route from 1T to 2T regimes is discontinuous, i.e., there is a jump, this is due to the presence of piecewise linear spring combined with time-varying stiffness, and it is typical of gear dynamics; this appears both for sun rotation, Fig. 7a, and first-planet rotation, Fig. 7b.
In the mesh frequency range ω m ∈ (5,750−7,000) Hz, case 2, Fig. 8 indicates the onset of gear tooth separation and chaotic response; because ω 4 = 2 × ω 1 , the third harmonic of ω 4 and the sixth harmonic of ω 1 , overlap with fourth harmonic of ω 2,3 = 1,445 Hz, and this is probably the reason of the chaotic region. Inside the interval ω m ∈ (7,000−7,200) Hz, close to a rotational mode resonance, a 2T region takes place for first-planet rotation although the motion is chaotic for sun rotation, see Fig. 8 and Table 9. Another chaotic region appears at ω m = ω 12 = 7,182 Hz. A sudden disappearance of a chaotic orbit is often called blue sky catastrophe; however, another motivation for a sudden change in dynamics can be the presence of coexisting stable chaotic and regular orbits, Ref. [37].
At frequency range ω m ∈ (13,400−13,800) Hz, Fig. 6, there is another chaotic region that leads to 4T region at the frequency range ω m ∈ (13,800−14,200) Hz. A further 2T region appears for ω m ∈ (14,200 − 14,800) Hz. Figures 9 and 10 give additional and extremely important information about the system dynamics for the case 2. Figure 10a corresponds to x-translation of sun gear center and Fig. 10b corresponds to its ytranslation. In linear field, the symmetry and the perfect balancing of the system implies that the sun is loaded with a self-equilibrate force system from the planets; therefore, it should experience no displacement in x and y directions. This is particularly evident in the case 2 Fig. 10. For the ranges where complex phenomena appear, the system presents symmetry breaking and consequently sun imbalance. Figure 11 presents the bifurcation diagram for all planets rotations. All planets undergo to the same vibration except for the chaotic regions when the symmetry breaking is evident.       Now the most interesting regimes found in the bifurcation diagrams are analyzed in detail. For each regime Poincaré maps, phase portraits, spectra and time histories are shown, see Table 9.
At 4,682 Hz (case 1), the dynamics appear complex, time histories appear nonstationary as well as the spectrum is dirty; however, the Poincaré map shows a scenario between quasiperiodic and strongly subharmonic; i.e., the distribution of points seems to lie on a continuous curve (typical of quasiperiodic), but the points are positioned in discrete locations. At 6,107 Hz (case 2), the gearbox experiences chaotic motion, the time histories and spectra are irregular; the phase trajectory fills a portion of phase space and the Poincaré map appears fractal.
At 6,808 Hz (case 2), there is probably a lowdimensional chaos, as a specific harmonic is prevalent in the spectrum.
As we have mentioned above and from Table 8, it is also clear at frequency range ω m ∈ (13,400−13,800) Hz that there is another chaotic region that leads to 4T region at the frequency range ω m ∈ (13,800−14,200) Hz. A further 2T region appears for ω m ∈ (14,200 − 14,800) Hz. It seems that here the route to the chaos is a cascade of period doubling (see Fig. 8).

Conclusions
This paper investigates the nonlinear dynamic behavior of a planetary gear system with equally spaced planets. The effect of mesh frequency is evaluated, and the system behavior at different rotational velocities is examined; this may help designers to understand how velocity can affect dynamic response of a plan-  Modes corresponding to multiples of mesh frequency can be excited due to: (1) the presence of higher harmonics of mesh stiffness variation and (2) nonlinear superharmonic resonance. Moreover, the resonance peaks lean to the left, implying softening nonlinearity induced by teeth separation.
Comparing the presents results with those of the work of Ambarisha and Parker, Ref. [2] shows that there are additional resonance peaks below 1,600 Hz, ∼ =1,800 Hz (the first distinct mode) and ∼ =8,000 Hz.
These peaks are the due to a combination of effects: (1) the parametric instability from higher harmonics of mesh stiffness and (2) nonlinear subharmonic and superharmonic resonances of the first and second distinct modes.
Reducing the translational bearing stiffnesses generally causes decreasing the natural frequency as we expect from elementary considerations. Moreover, changing the bearing stiffness can influence the ratios between modes and, for certain values, produce additional internal resonances; this is reflected in a higher modal interaction in nonlinear regimes.
Different nonlinear phenomena such as nonlinear jumps, chaotic motions and period-doubling bifurcations occur when the mesh frequency or any of its higher harmonics is near a natural frequency of the system. The nonlinearity is due to teeth separation that take place when the amplitude of vibration generates inertia forces exceeding the static preload.
The occurrence of a parametric instability is found when the rotational speed is twice the resonance of a rotational mode; this phenomenon is typical of almost all gearboxes.
In the presence of internal resonances, remarkable changes are observed; in particular, the variation of the bearing stiffness can greatly influence and enrich the dynamic scenario.
One of the most interesting disclosures of this work is related to the chaos-induced imbalance; i.e., due to the symmetry of the system, there is a perfect force balance; this means that the sun bearings are subjected to negligible forces (theoretically zero), and this is true in principle both in linear and nonlinear regular vibrations. However, when the vibration is chaotic, the system experiences a symmetry breaking, a wellknown phenomenon in chaotic dynamics. This symmetry breaking in the response has an important counterpart in the mechanical behavior of the system, i.e., it generates unexpected and undesirable loads on the sun supports. This aspect of the system dynamics cannot be accounted for using the classical design tools, leading to possibility of new failure modes.