Bifurcation Observation of Combining Spiral Gear Transmission Based on Parameter Domain Structure Analysis

This study considers the bifurcation evolutions for a combining spiral gear transmission through parameter domain structure analysis. The system nonlinear vibration equations are created with piecewise backlash and general errors. Gill’s numerical integration algorithm is implemented in calculating the vibration equation sets. Based on cell-mapping method (CMM), twodimensional dynamic domain planes have been developed and primarily focused on the parameters of backlash, transmission error, mesh frequency and damping ratio, and so forth. Solution demonstrates that Period-doubling bifurcation happens as the mesh frequency increases; moreover nonlinear discontinuous jump breaks the periodic orbit and also turns the periodic state into chaos suddenly. In transmission error planes, three cell groups which are Period-1, Period-4, and Chaos have been observed, and the boundary cells are the sensitive areas to dynamic response. Considering the parameter planes which consist of damping ratio associatedwith backlash, transmission error,mesh stiffness, and external load, the solution domain structure reveals that the system step into chaos undergoes Period-doubling cascade with Period-2 (m: integer) periodic regions. Direct simulations to obtain the bifurcation diagram and largest Lyapunov exponent (LE) match satisfactorily with the parameter domain solutions.


Introduction
The geared transmission which is driven by several powers is considered as a combining geared system.In contrast to the ordinary gear structures, combining configuration can readily transmit the power from separate branches into one route to deliver a strongly transmissible drive.Their inherent distinctive features make it fit very well in some specific fields.One of its significant applications is in the cabin of today's mainstream helicopters, such as Blackhawk UH-60, Apache AH-64, Tiger EC-665, and Mil Mi-26, as the requirement of high performance working condition needs deep-going investigation on the gear dynamic characteristics.Generally, the gear system is subjected to various impulsive excitations including backlash, geometry error, time-varying mesh stiffness, and torque fluctuations.These parametric excitations exist in the driving situation and relate to the tooth engagement stability and some significant nonlinear instabilities such as bifurcation and chaos.For the nonlinear system, bifurcation is a typical phenomenon and has been conducted through various techniques.Poincaré map, phase portrait, and Largest Lyapunov exponent already have been exhibited efficiently in dealing with the matters with bifurcation parameter, but when the system is concerned with multiple quantities, the long term prediction of the global analysis to the system has to take other strategies to explore.The observation of the global bifurcation is necessarily relying on the solutions of the dynamical domains.In the parameter plane, the dynamic domain structure provides a visible sight to the routes opening into chaos, and the size of dynamic group reveals the effect of the control parameters.According to the cell state, the bifurcation solutions are a combination of multiple dynamic motions; once the bifurcation parameter varies, the system will respond with a change of the corresponding attractor.Therefore, a better way to investigate the bifurcation behavior is to trace the development of the dynamical domains.Since the cell-mapping method has been successfully applied in the field of global exploration, here it is considered to divide the parameter plane into series of excitation cells to locate the existence of all the attractors, and based on these one can determine the bifurcation developments.
As for the gear system, geometry clearance, transmission error, and variation meshing stiffness are of interest.Such vibratory elements may interact mutually to bring about the unexpected dynamical phenomena.In the early literature [1][2][3][4], the appearance of the saddle-node bifurcation, grazing bifurcation, and Hopf bifurcation leading into chaos has been presented.Researchers have discovered homoclinic bifurcation and intermittent chaos under the excitations of instability domains as well [5][6][7][8][9][10][11].Lin and Parker [12] examined parametric instability in view of fluctuating stiffness under multiple mesh conditions.Li et al. [13] systematically investigated bifurcation and chaos properties in terms of damping ratio, backlash, and excitation frequency and found the system motion state changes into chaos via Hopf bifurcation finally.Studies by Ambarisha and Parker [14] predicted that, due to the unstable contact loss and periodically timevariable characters, bifurcation and chaotic perturbation will occur when the mesh frequency is close to a natural frequency.It is obvious that smooth running, quiet noise, and long life service are an expectation in the gear actual applications.However, the unpredicted problems such as displacement shock and intermittent perturbation may happen on the meshing tooth when the parameter changes.In order to identify the dominant excitation related to the existing issues, many researchers have been dedicated to the explorations.Chang-Jian and Chang [15] examined the chaotic behavior with respect to the sensitivity of participation parameters; the systematic investigation was studied in conjunction with Lyapunov exponent, fractal dimension, Poincaré section, and bifurcation diagrams.With a single degree of freedom geared model, Sato et al. [16] have considered the harmonic excitations and investigated the chaotically transitional phenomena through Li-Yorke's theorem, and the tangent bifurcation sets and pitchfork bifurcation sets were plotted on the parametric plane.Parker and Guo [17] based on the planetary gear nonlinearity further analyzed the solution stability using Floquet theory.Their simulation reveals that the secondary Hopf bifurcation was referred to the transition from quasi-chaos to chaos.
In the prior literatures, one can find that bifurcation characteristics and nonlinear behaviors have been discussed through several approaches.Nevertheless, a deeper understanding of the bifurcation evolution as well as the transition between stable periodic and chaotic motion is still essential for designing and adjusting the key parameters to enable the system to perform well.Since Hsu first proposed the cellmapping method to cope with the global analysis for the nonlinear system, it has proven to be an efficient technique to locate the strange attractors as well as basins of attraction [18,19].With the development of dynamical technologies, more and more attentions on the global analyses have been given.Hinrichs et al. [20] examined a nonsmooth oscillator.The solution about the domain of attraction has given insight into the overall dynamic behavior; meanwhile, the special attention to the discontinuity region was also pointed out.de Souza et al. [21] have focused on the researches of parameter intervals for a rattling gearbox system and found that the basin hopping and different attractors switching are due to the additional noise.Yunwen et al. [22] developed a mixed cell mapping method; the obtained coexisting attractors suggested that parameter structure has the significant advantage in performing gear dynamic behaviors.
Most recently, Farshidianfar and Saghafi employed Melnikov analytical approach to trace the global homoclinic bifurcation along with the evolution of chaotic orbit.Their analytical simulations demonstrate the homoclinic bifurcation and the transition route to chaos with respect to the stable and unstable manifolds [23].The remainder of this paper is organized as follows.In Section 2 dynamic model involving nonlinear backlash is formulated and solved utilizing step-variant Gill numerical algorithm.In Section 3.3 the bifurcation evolution procedure is conducted and observed with respect to transmission error and mesh frequency.In Section 3.4 two-dimensional parameter domain structure is explored, and the bifurcation diagram and largest Lyapunov exponent are used to validate the solutions.Finally, Section 4 draws conclusions and describes the outcome of the simulation results.

Dynamics Model
In Figure 1, for a combining spiral gear system, the shaft angle between mating gears is set as /2.We introduce the generalized coordinate for gear modeling, where the origins   stand for the geometric centre of each gear.Under the assumption of a rigid gear body, the dynamical model can be constructed by employing a lumped-parameter method, where the gear body is described by mass   , and with equal support bearing stiffness and damping in , ,  directions.In the normal direction of tooth profile, the backlash and static transmission error are considered.Since each deflection displacement can be described by means of the degrees of freedom (DOF), one component mass would possess 3 translational DOF, namely, , , and , and 1 angular rotation, which can be described by   =     , where   and   are the meshing point radius and angular displacement; eventually, the system would amount to be of 12 DOF totally.Throughout In order to eliminate the rigid body displacement, the normal relative displacements  13 and  23 have to be investigated by means of the general coordinates.Thereupon we assume that the deflection along the line of action is positive; accordingly, the translational component deflection   ,   , and   and angular deflection   can be formulated by Δ 1 , Δ 2 , Δ 3  1 , and Δ 3 2 .Similarly, let the mesh forces  13 and  23 project on the coordinate axis with three distributed forces, which can be expressed by   ,   , and   .In Table 1, the main parameters of the study gear system are listed.
Due to the symmetry structure of pinion 1 and pinion 2, the axial component displacements and component forces projected in the meshing direction can be given by where Δ 1 and Δ 2 designate the relative displacement vector for pinion 1 and pinion 2, respectively, and with respect to their force vectors  1 and  2 , respectively.
A similar operation can be done on driven gear 3, given by where Δ 3 1 and Δ 3 2 represent the relative displacement vectors for gear 3 and  3  1 and  3 2 denote the force vectors acting on gear 3, where the subscripts 1 and 2 denote the left side and right side, respectively.Now, it is apparently that the overall component force acting on gear 3 becomes In ( 1)-( 3), the term   ( = 1, 2, . . ., 5) is used repeatedly in calculating the mesh force and the deflection displacement is a derivative parameter.Expressions are given as follows: The translational equations of motion for the system can be obtained as follows: where   = {  ,   ,   }   = 1, 2, 3 stand for the translational degrees of freedom.
The differential equations for torsional motion for each body in the direction of   are formulated as follows: Equations ( 5)-( 6) are the system governing equations of motion and can be modeled as a set of second-order differential equations: where The meshing force   acting between contact teeth is where  = 1, 2,  = 3,    () =    (1 +  sin(Ω +   )),    is the mean mesh stiffness between gear pair,  is the stiffness factor, and the mesh frequency and initial mesh phasing are represented by Ω and   , respectively.The relative deflection   is coupled with the translational and torsional motions, associated with the transmission error, and can be expressed as where  = 1, 2,  = 3;   () is the time varying general error, expressed by   () =  sin(Ω +   ), and  is the amplitude of transmission error.
Substituting the quantities of Δ 1 , Δ 2 , Δ 3 1 , and Δ 3 2 into (9), the coupled formulation can be rewritten as follows: where Eliminating rigid body motion (7) then can be rewritten as where q(t) = { 1 ,  Other terms in (12) are listed below: Mathematical Problems in Engineering 5 where and   are the equivalent masses introduced to simplify the computational process, defined by Mesh damping    associated with the mesh stiffness can be obtained by where   is the damping ratio (  = 0.01).
Nonlinear backlash function is described by a piecewise formulation (,   ), given by To enhance the ability of convergence in the differential equation computation, we introduce the terms    and   for the nondimensionalization operation, where    = 100 × 10 −6 m represents a characteristic length and   is the characteristic frequency used to define dimensionless time  as well as other quantities.Here   is defined by With these terms, (12) yields the dimensionless equations of motion and is rewritten in the following forms: where ζ = {  ,   ,   }  ,  = 1, 2, 3.

Bifurcation Behaviors
For ( 12), as a nonautonomous equation set with multiple parameters, it can be generally described as follows: where (, , ) is a Jacobi matrix,  stands for the excitation source in the system, and  0 represents an initial condition.

Lyapunov Exponent. Introduce two close initial values
5 ( 13 (0), δ 13 (0)) and   5 (  13 (0), δ  13 (0)), and define them by  0 = ‖ 5 −    5 ‖ = 0.0001, where  0 stands for the distance of those two state vectors.After several steps of computation, the initial value would become  5 ( 13 (), δ 13 ()) and    5 (  13 (), δ  13 ()), and the separation of those two orbits can be evaluated through [24,25] While (23) reaches the steady state, the state variable would develop into  5 () and   5 () at the time step of ()th, where  is the time step in Gill's computation,  is the current number of iteration, and the distance at this moment is With the iteration of the state points, the result on the orbits at the next step ( + 1)th would be  5 [( + 1) ] =    5 () , where   represent the mapping formula under parametric excitation.Combining ( 24) and ( 25) the trajectory separation at ( + 1)th step will be One can find a start point   5 that will keep a constant distance  0 away from  5 [( + 1)] which avoids the divergence of  +1 during numerical procedure.After a certain number of steps, a series of Lyapunov exponents are obtained as follows: where  = 1, 2, . . ., ∞.
The largest Lyapunov exponent  1 in terms of ( 27) can be obtained by 3.2.Cell-Mapping Method.In order to reveal the global bifurcation behavior in parametric space, we adopt the CMM [26][27][28][29][30].In this study the global domain comprised several control quantities from the system equations and is uniformly described by   = {  }  , where   is one of the concerned parameters, in the global domain with the forms of where   is a unit cell in domain space.
The centre of the cell is used in computing.For an dimensional domain, the cell identification can be discretized as where  refers to the quantities of interest and  refers to the th cell.
Here, only a two-dimensional domain is used to simplify the calculation for convenience, which also satisfies the requirement of gear dynamic exploration.In this case, the centre of the th cell is transformed as Under the periodical excitation, the effect with regard to   has the following relation to the system response: Data is only captured after the computation runs into steady state, according to (32), and supposing the system is in the Period-1 state, it can be described by  +1 = (  ,   ), if it is a multiple periodic motion, and the mapping route of term   would be (  ,   ) →  (  ,   ) →  2 (  ,   ) → ⋅ ⋅ ⋅ →   (  ,   )  ∈ . (33) If  → ∞, (33) is still not satisfied, which means that the cell of   is in a chaotic state.
In a dynamic region Σ which contains  ×  cells, the whole domain may consist of the periodic, quasi-periodic, or chaotic state due to the analytic solution.If the cells   demonstrate being -periodic motion, such cells can be described by The union of quasi-periodic cell is Similarly, it yields the chaotic domain The whole domain plane can be assembled by

Bifurcation Evolution.
With respect to the bifurcation investigation, the key parameter   is set in a certain range and varies with time ; accordingly a series of state vectors can be obtained via   (  , ),  = 1, 2, . . ., .Here, the bifurcation diagram is considered using deflection displacement  13 versus variable   which can be plotted by (  ,  13 ).Figure 2 demonstrates the threshold of bifurcation as well as the development process from a simple motion to chaos with the changes of .At first  = 28, three bifurcation nodes present at  1 = 0.917,  2 = 0.935, and  3 = 0.951; within  ∈ [0.9, 0.954], the system is completely under periodic state without uncertainty.It is also shown that a discontinuous jump is taking place at the second node, and this break makes the Period-2 orbit be split as two sections before coming into Period-1 motion; meanwhile, the instability region emerges from  > 0.954.When  = 28.and more degenerate and accompanied by further and further strengthening of the unstable vibration at  > 0.951.For  = 30 and  = 32, the right side of the Period-2 orbit has disappeared, giving way to stable Period-1 vibration instead.
Along with the bifurcation development, the nonlinearity manners become stronger gradually under the combination effects of error and mesh frequency, which makes the system response in turn alternate between the periodic and chaotic motions.At  = 42, we see the initial right branches of Period-2 which have been totally replaced by Period-1 motion over the span of  ∈ [0.9, 0.96].Finally, in Figure 2(f) the bifurcation diagram presents three evident periodic windows at  1 ∈ [0.975, 0.978],  2 ∈ [0.989, 0.999], and  3 ∈ [1.012, 1.016], which have been verified with periodic responses of Period-9, Period-11, and Period-11, respectively.Moreover, the amplitude of bifurcation diagram is increased further too.One thing to be noted is that the occurrence of the bifurcation intersection location shifts to a larger value of the mesh frequency in the process of error increasing.This may be considered as an opportunity to improve the system bifurcation characteristics.Displacement bifurcation versus dimensionless frequency has been shown in Figure 3 One can find two significant skips that appear at  = 1.927 and  = 2.112.We note that the first one also breaks the route of Period-2 orbit, and the other one makes the system enter chaos with a suddenly induced exterior crisis due to the collision with a chaotic attractor, both of which are more likely to induce instability or tooth collision in gear system.As expected, the largest Lyapunov exponent in Figure 3(b) provides a good determination regarding the transition from periodicity to chaos in the light of the criterion value of 0.

Domain Structure Analysis.
Discretize the parameter space to be a two-dimensional plane with   ⊗   , where   denotes the th cell.The mapping initial condition is set as  5 = 0, letting  1 = 1 − 2 and  2 = 1 − 4 to be the computing precision to distinguish the target cell unit and identify the dynamic periodicity; for instance, if | 5 ()− 5 (+1)| <  1 , that means that the dynamic state is under periodic motion, then judging the periodicity by comparing with  2 .In order to capture the accurate solutions, we record 1000 mapping periods for a single cell parameter.
In Figure 4(a) the whole plane is divided into 22500 rectangular cells with uniform size of intervals of 0.008 and 0.0008 for  and   , respectively.Clearly, the global attraction is constructed with three distinct domains of 1, 4, and Chaos, where  represents the periodic domains; the numbers 1, 4 refer to the motion characteristics.It is obvious that the window of Period-4 orbit is filled inside the chaotic domain.With increasing the transmission error , the periodic motion and chaotic state present alternately in the parameter plane.In Figure 4 group mainly consists of 1, 2, and 3 as well as Chaos, with some regular cells forming an instability narrow zone, which indicates that the system behavior is sensitive with the impact of errors and mesh frequency.It could be convinced that this area is unstable under the corresponding excitation, for the gears should avoid being in such irregular region.In this ⊗ plane, Period-2 cells have a spatial domination over the system response.During the changes of decreasing , the Period-1 cells disappeared firstly and the Period-3 orbit moves upwards along the  edge on the left side of Period-2 group.Moreover, over the range of  ∈ [0.74, 0.825], the cells mapping into chaos is not smooth.As  goes up to 0.82, the cell evolution undergoes a transition from Period-1 motion to an irregular region and then to chaotic domain, and chaos appears following Period-3 cells.
In Figure 6(a), the occurrence of bifurcation motion is synchronous with largest Lyapunov exponent in the case  = 0.788, and the system vibration goes through stable Period-1 and Period-2 orbit to the chaotic vibration.Bifurcation diagram and the largest LE both discover that a sudden jumping occurs at  = 1.45 just happening at Period-2 motion which is transiting into chaos.At this moment, intermittent fluctuation also can be found from the largest Lyapunov exponent with a sharp switch from −0.02 to 0.02.This unstable parameter domain probably shadows an existence of chaotic domain nearby.However, CMM does not have the capability to find such a nonlinear jumping discontinuity; yet it gives a satisfying conformation about the emergence of chaos.
For Figure 7, the CMM solutions are depicted employing the damping ratio of  ∈ [0.05, 0.41] and, respectively, versus parameters of backlash   , mesh error , stiffness factor , and input torque  1 .In these four plots, likewise, chaotic cells cover a big area with small  inside the parametric planes.Period-doubling bifurcation is exhibited along the route approaching chaos and passes in sequence through 2, 4, 8, and 16; besides, some 32 points are also evident at the edge before chaos emerges.In Figures 7(a)-7(c) after Period-4 mapping region, small quantity of Period-8 and Period-16 cells can be seen.In  ⊗  1 plane, very few Period-32 points are showing up at the intersection near the onset of chaos.The reason is that, with the bifurcation of the periodic motions, the windows of the coming periodic response become narrower and narrower, and small cells are able to explore the transition.In Figure 7(a), three Pseudoperiodic cells are noticed lying on a deviate location, which fails to reveal the real nature of the dynamic motions.Further investigation of the periodic evolvement is subject to a finer resolution of parametric field.According to these four planes, we know that the mapping result shows specific structure in the evolution plane which indicates various bifurcation manifestations for the system.Specifically, the domains of attraction gather to develop special features as parameter changes, which substantially helps to track different dynamic responses.It is also interesting that the system is in the chaotic state for a small  and in periodic motion at larger .Hence one should take into account that a higher damping ratio will suppress the gear vibration.
The condition in Figure 8 is pictured with 300 steady periods for a fixed .With the damping ratio decreasing, it mainly illustrates the changes starting from a Period-2  response (: integer) to a doubled periodic ( + 1) again and goes on until eventually reaching chaos with the recursive topological structure, and we can find that the mesh movement experienced a dramatic transition.Figure 8(b) shows that the largest Lyapunov exponent rises gradually from the value of −0.1 to 0.08 ultimately, where the largest LE indeed matches the tendency of bifurcation development.

Conclusions
The nonlinear vibration equations of combining spiral gear transmission were set up, and the calculations were done by adopting Gill's numerical algorithm.The investigation considered gear backlash, transmission error, damping ratio, and variation mesh stiffness.Using the cell-mapping technique leads to good insight into the dynamic features of parameter domain structures.Bifurcation diagram and Lyapunov exponent are used to validate the results in parameter spaces.Some main conclusions are as follows: (1) Under the excitation of mesh frequency , perioddoubling bifurcation was observed with the Largest Lyapunov exponent  1 ∈ [−0.021, 0.038]; when  = 1.93, the discontinuous jump breaks the Period-2 orbit and at the point of  = 2.12, the jump switches Period-2 state into chaos suddenly.As the transmission error changes from 28 to 42, bifurcation evolution was exhibited.
(2) In  ⊗   and  ⊗  parameter spaces, three dynamic regions including Period-1, Period-4, and Chaos were exhibited, and gear dynamic responses are sensitive to the boundary parameters of the cell groups.
(3) In  ⊗   ,  ⊗ ,  ⊗ , and  ⊗  1 planes, the parameter domain structure reveals the global dynamic behavior distributions, with periodic regions of -2, -4, -8, . .., and Period-2  (: integer), leading to chaos, and large damping ratio contributes to a steady dynamic response.Initialmeshphasing Ω: Mesh frequency with dimension : Dimensionless amplitude of errors   : Th ee q u i v a l e n tm a s s e s   : Th ee q u i v a l e n tm a s so fg e a r s  =   /  : Th em e s hd a m p i n gr a t i o : Dimensionless mesh frequency : Displacement vector M: Total mass matrix K: Total stiffness matrix C: Total stiffness matrix a  : Computation coefficient sub-matrix.

3 Figure 1 :
Figure 1: Dynamic model for a combining spiral gear transmission using lumped-parameter method.

Figure 6 :
Figure 6: Verification of the domain structure with respect to transmission error.

Figure 7 :
Figure 7: Domain structure analysis in two-dimensional parameter planes.

Figure 8 :
Figure 8: Bifurcation evolution with respect to the damping ratio.

Table 1 :
Main dynamic parameters for the study gear system.
1 ,  1 ⏟⏟⏟⏟⏟⏟⏟⏟⏟⏟⏟⏟⏟⏟⏟ ,  13 ,  23 }  , q(t) is the new degrees of freedom of the system, and  13 and  23 are no longer containing rigid body displacement.isthehalf backlash.Obviously, (12) turns out to have 11 DOF after simplification.