#### Abstract

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), two-dimensional 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 associated with 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*- (: integer) periodic regions. Direct simulations to obtain the bifurcation diagram and largest Lyapunov exponent (LE) match satisfactorily with the parameter domain solutions.

#### 1. 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–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–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 time-variable 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 cell-mapping 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.

#### 2. Dynamics Model

In Figure 1, for a combining spiral gear system, the shaft angle between mating gears is set as . 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 this paper, if there is no special designation, the subscript denotes pinion 1, pinion 2, and bevel 3, respectively, and superscript represents the coordinate axis of , , and .

In order to eliminate the rigid body displacement, the normal relative displacements and 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 , , , and . Similarly, let the mesh forces and 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 bywhere and designate the relative displacement vector for pinion 1 and pinion 2, respectively, and with respect to their force vectors and , respectively.

A similar operation can be done on driven gear 3, given bywhere and represent the relative displacement vectors for gear 3 and and denote the force vectors acting on gear 3, where the subscripts and 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 () 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 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 iswhere , , , 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 aswhere , ; is the time varying general error, expressed by , and is the amplitude of transmission error.

Substituting the quantities of , , , and into (9), the coupled formulation can be rewritten as follows:where

Eliminating rigid body motion (7) then can be rewritten aswhere , is the new degrees of freedom of the system, and and are no longer containing rigid body displacement. is the half backlash.

Obviously, (12) turns out to have 11 DOF after simplification.

Other terms in (12) are listed below:where

and are the equivalent masses introduced to simplify the computational process, defined by

Mesh damping associated with the mesh stiffness can be obtained bywhere is the damping ratio ().

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 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 , .

Other essential nondimensional parameters used are defined as follows: means dimensionless unit.

The varying-step Gill integrations are applied to compute the dimensionless equation sets. For the convenience of programming, define the following state quantities, , , , , and **,** to make the system equations become a first-order differential equation set, given bywhere , , represent the translational displacement vector and () represent the relative displacements.

#### 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 represents an initial condition.

##### 3.1. Lyapunov Exponent

Introduce two close initial values and , and define them by , where stands for the distance of those two state vectors. After several steps of computation, the initial value would become and , 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 and 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 th would bewhere represent the mapping formula under parametric excitation. Combining (24) and (25) the trajectory separation at th step will be

One can find a start point that will keep a constant distance away from which avoids the divergence of during numerical procedure. After a certain number of steps, a series of Lyapunov exponents are obtained as follows:where .

The largest Lyapunov exponent 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–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 ofwhere 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 aswhere 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 , if it is a multiple periodic motion, and the mapping route of term would be

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

##### 3.3. 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 , . Here, the bifurcation diagram is considered using deflection displacement versus variable which can be plotted by . 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 , three bifurcation nodes present at , , and ; within , 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 . When and , the* Period-*2 motion over the range becomes more and more degenerate and accompanied by further and further strengthening of the unstable vibration at . For and , 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 , we see the initial right branches of* Period-*2 which have been totally replaced by* Period-*1 motion over the span of . Finally, in Figure 2(f) the bifurcation diagram presents three evident periodic windows at , , and , 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.

**(a)**

**(b)**

**(c)**

**(d)**

**(e)**

**(f)**Displacement bifurcation versus dimensionless frequency has been shown in Figure 3(a) based on . One can find two significant skips that appear at and . 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.

**(a) Bifurcation diagram**

**(b) Largest Lyapunov exponent**

##### 3.4. 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 , letting and to be the computing precision to distinguish the target cell unit and identify the dynamic periodicity; for instance, if , that means that the dynamic state is under periodic motion, then judging the periodicity by comparing with . 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 , , and* Chaos*, where represents the periodic domains; the numbers 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(c) the control parameters for the bifurcation diagram are extracted from the plane in the case and , and we can see that the response suddenly alters at the bifurcation point, which means a slight disturbance close to the intersection of this parametric area would result in significant breaking of the equilibrium state, with original attractor destroyed suddenly. The cells in show that backlash plays a small role on the bifurcation adjustment when . Another parameter domain plane with was conducted in Figure 4(b) and received the similar attractions with , , and* Chaos*. Bifurcation diagram in Figure 4(d) indicates a good agreement with the parameter domain analysis.

**(a)**

**(b)**

**(c)**,

**(d)**,The results in Figure 5 demonstrate the complicated bifurcation evolutions in the bounded domain. The dynamic group mainly consists of , , and 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 , the cells mapping into chaos is not smooth. As goes up to , 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 , 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 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.

**(a)**Bifurcation diagram with**(b) Largest Lyapunov exponent**

For Figure 7, the CMM solutions are depicted employing the damping ratio of and, respectively, versus parameters of backlash , mesh error , stiffness factor , and input torque . 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 , , , and ; besides, some 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 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 Pseudo-periodic 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.

**(a)**

**(b)**

**(c)**

**(d)**The condition in Figure 8 is pictured with steady periods for a fixed . With the damping ratio decreasing, it mainly illustrates the changes starting from a* Period*- response (: integer) to a doubled periodic () 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.

**(a)**Bifurcation diagram with**(b) Largest Lyapunov exponent**

#### 4. 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 ,* period*-doubling bifurcation was observed with the Largest Lyapunov exponent ; when , the discontinuous jump breaks the* Period*-2 orbit and at the point of , 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 planes, the parameter domain structure reveals the global dynamic behavior distributions, with periodic regions of , and* Period*- (: integer), leading to chaos, and large damping ratio contributes to a steady dynamic response.

#### Nomenclature

: | Mass of the gear- |

: | Gear transmission ratio |

: | Support bearing stiffness |

: | Support bearing damping |

, : | Gear meshing stiffness |

, : | Gear meshing damping |

, : | Input torques on pinions |

: | Output torque on gear 3, where |

: | Moment of inertia for each gear |

, : | Relative displacement functions |

: | Generalized coordinate axis |

: | Origin of coordinate system |

: | Reference cone angle |

: | The half backlash |

, : | Displacement vector for pinions |

, : | Vibration displacement vector on gear 3 |

, : | Mesh force acting at tooth contact point |

: | Component forces in axis direction |

: | Computation coefficient |

, : | Gear transmission error functions |

, : | Initial mesh phasing |

: | Mesh frequency with dimension |

: | Dimensionless amplitude of errors |

: | The equivalent masses |

: | The equivalent mass of gears |

: | The mesh damping ratio |

: | Dimensionless mesh frequency |

: | Displacement vector |

: | Total mass matrix |

: | Total stiffness matrix |

: | Total stiffness matrix |

: | Computation coefficient sub-matrix. |

*Superscript*

: | . |

*Subscript*

: | . |

#### Competing Interests

The authors declare that there is no conflict of interests regarding the publication of this paper.

#### Acknowledgments

This study is based upon financial supported by the National High Technology Research and Development Program of China (2009AA04Z404). The first author also acknowledges the grateful discussion with Professor Earl H. Dowell at Duke University.