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

Authors: He Lin; Sanmin Wang; Earl H. Dowell; Jincheng Dong
Journal: Mathematical Problems in Engineering
 (2016)

---

## 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 thatPeriod-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 arePeriod-1,Period-4, andChaos 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 undergoesPeriod-doubling cascade withPeriod-2
m (m: integer) periodic regions. Direct simulations to obtain the bifurcation diagram and largest Lyapunov exponent (LE) match satisfactorily with the parameter domain solutions.


---

## Body

## 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 Section2 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 Figure1, 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 o
i 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 m
i, and with equal support bearing stiffness and damping in x, y, z 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, x, y, and z, and 1 angular rotation, which can be described by u
i
=
r
i
θ
i, where r
i and θ
i 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 i
=
1,2
,
3 denotes pinion 1, pinion 2, and bevel 3, respectively, and superscript j represents the coordinate axis of x, y, and z.Figure 1
Dynamic model for a combining spiral gear transmission using lumped-parameter method.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 x
i, y
i, and z
i and angular deflection θ
i can be formulated by Δ
1, Δ
2, Δ
1
3, and Δ
2
3. Similarly, let the mesh forces F
13 and F
23 project on the coordinate axis with three distributed forces, which can be expressed by F
x, F
y, and F
z. In Table 1, the main parameters of the study gear system are listed.Table 1
Main dynamic parameters for the study gear system.
Items
Pinion 1 (2)
Gear 3
Number of teethz
23
67
Modulem (mm)
3
3
Pressure angleα
n (deg)
25
25
Helix angleβ
i (deg)
35
35
Massm
i (kg)
2.6
6.8
Moment of inertiaJ
i (kg⋅m2)
0.018
0.042
Bearing stiffnessk
i
j (N/m)
2.08 × 108
2.08 × 108
Mean mesh stiffnessk
m
1, k
m
2 (N/m)
1.21 × 109
1.21 × 109Due 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(1)Δ1=-a2x1-a3y1-a1z1a3r1θ1,Δ2=a3x2a2y2-a1z2a3r2θ2,F1=a2F13a3F13a1F13,F2=-a3F23-a2F23a1F23,where Δ
1 and Δ
2 designate the relative displacement vector for pinion 1 and pinion 2, respectively, and with respect to their force vectors F
1 and F
2, respectively.A similar operation can be done on driven gear 3, given by(2)Δ13=a3x3-a5y3a4z3-a3r3θ3,Δ23=-a3x3a5y3a4z3-a3r3θ3,F13=-F13a3F13a5-F13a4,F23=F23a3-F23a5-F23a4,where Δ
1
3 and Δ
2
3 represent the relative displacement vectors for gear 3 and F
1
3 and F
2
3 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(3)F3=-F13a3F13a5-F13a4+F23a3-F23a5-F23a4.In (1)–(3), the term a
i (i
=
1,2
,
…
,
5) is used repeatedly in calculating the mesh force and the deflection displacement is a derivative parameter. Expressions are given as follows: (4)a1=sin⁡γ1sin⁡αn+cos⁡γ1sin⁡β1cos⁡αn,a2=cos⁡γ1sin⁡αn-sin⁡γ1sin⁡β1cos⁡αn,a3=cos⁡β1cos⁡αn,a4=cos⁡γ3sin⁡β3cos⁡αn-sin⁡γ3sin⁡αn,a5=cos⁡γ3sin⁡αn+sin⁡γ3sin⁡β3cos⁡αn.The translational equations of motion for the system can be obtained as follows:(5)m1∂2ζ1∂t2+c1ζ1∂ζ1∂t+k1ζ1ζ1+a1ζ1cm1δ˙13+km1tfδ13=0,m2∂2ζ2∂t2+c2ζ2∂ζ2∂t+k2ζ2ζ2+a2ζ2cm2δ˙23+km2tfδ23=0,m3∂2ζ3∂t2+c3ζ3∂ζ3∂t+k3ζ3ζ3+a3ζ3cm1δ˙13+km1tfδ13+a3′ζ3cm2δ˙23+km2tfδ23=0,where ζ
i
=
x
i
,
y
i
,
z
i
T
i
=
1,2
,
3 stand for the translational degrees of freedom.The differential equations for torsional motion for each body in the direction ofθ
i are formulated as follows:(6)J1∂2θ1∂t2=T1-Fx1r1,J2∂2θ2∂t2=T2+Fx2r2,J3∂2θ3∂t2=-T3-Fx3r3+Fx3r3.Equations (5)-(6) are the system governing equations of motion and can be modeled as a set of second-order differential equations:(7)Mq¨+Cq˙+Kq=Fij,where {
q
}
=
x
1
,
y
1
,
z
1
,
θ
1
︸
pinion
1
,
x
2
,
y
2
,
z
2
,
θ
2
︸
pinion
2
,
x
3
,
y
3
,
z
3
,
θ
3
︸
bevel
3
T.The meshing forceF
i
j acting between contact teeth is(8)Fij=kmitfδij+cmiδ˙ij,where i
=
1,2, j
=
3, k
m
i
(
t
)
=
k
m
i
(
1
+
m
k
sin
⁡
(
Ω
t
+
ψ
i
)
), k
m
i is the mean mesh stiffness between gear pair, m
k is the stiffness factor, and the mesh frequency and initial mesh phasing are represented by Ω and ψ
i, respectively. The relative deflection δ
i
j is coupled with the translational and torsional motions, associated with the transmission error, and can be expressed as(9)δij=∑ifxi,yi,zi,ui+∑jfxj,yj,zj,uj-eit,where i
=
1,2, j
=
3; e
i
(
t
) is the time varying general error, expressed by e
i
t
=
E
n
sin
⁡
(
Ω
t
+
ψ
i
), and E
n is the amplitude of transmission error.Substituting the quantities ofΔ
1, Δ
2, Δ
1
3, and Δ
2
3 into (9), the coupled formulation can be rewritten as follows:(10)δ¨13+a2x¨1+a3y¨1+a1z¨1-a3r4∂2θ1∂t2-a3x¨3+a5y¨3-a4z¨3+a3r3∂2θ3∂t2=-e¨1t,δ¨23-a3x¨2-a2y¨2+a1z¨2+a3r2∂2θ2∂t2+a3x¨3-a5y¨3-a4z¨3-a3r3∂2θ3∂t2=-e¨2t,where(11)∂2θ1∂t2=T0-a3r1F13J1,∂2θ2∂t2=T0-a3r2F23J2,∂2θ3∂t2=-T1+a3r3F13+a3r3F23J3.Eliminating rigid body motion (7) then can be rewritten as(12)Mq¨t+Cq˙t+Kqft,bn=Ft,where q
(
t
)
=
x
1
,
y
1
,
z
1
︸
poinion
1
,
x
2
,
y
2
,
z
2
︸
pinion
2
,
x
3
,
y
3
,
z
3
︸
bevel
3
,
δ
13
,
δ
23
T, q
(
t
) is the new degrees of freedom of the system, and δ
13 and δ
23 are no longer containing rigid body displacement. b
n is the half backlash.Obviously, (12) turns out to have 11 DOF after simplification.Other terms in (12) are listed below:(13)
M
=
m
1
0
3
×
3
0
3
×
3
0
1
×
1
0
1
×
1
0
3
×
3
m
2
0
3
×
3
0
1
×
1
0
1
×
1
0
3
×
3
0
3
×
3
m
3
0
1
×
1
0
1
×
1
α
1
0
1
×
3
α
2
1
0
1
×
1
0
1
×
3
a
3
α
4
0
1
×
1
1
,
K
=
k
1
ζ
0
3
×
3
0
3
×
3
a
1
ζ
k
m
1
t
0
3
×
1
0
3
×
3
k
2
ζ
0
3
×
3
0
3
×
1
a
2
ζ
k
m
2
t
0
3
×
3
0
3
×
3
k
3
ζ
a
3
ζ
k
m
1
t
a
3
′
ζ
k
m
2
t
0
1
×
3
0
1
×
3
0
1
×
3
a
3
2
k
m
1
t
M
e
a
3
2
k
m
2
t
M
3
0
1
×
3
0
1
×
3
0
1
×
3
a
3
2
k
m
1
t
M
3
a
3
2
k
m
2
t
M
e
,
C
=
c
1
ζ
0
3
×
3
0
3
×
3
a
1
ζ
c
m
1
0
3
×
1
0
3
×
3
c
2
ζ
0
3
×
3
0
3
×
1
a
2
ζ
c
m
2
0
3
×
3
0
3
×
3
c
3
ζ
a
3
ζ
c
m
1
a
3
′
ζ
c
m
2
0
1
×
3
0
1
×
3
0
1
×
3
a
3
2
c
m
1
M
e
a
3
2
c
m
2
M
3
0
1
×
3
0
1
×
3
0
1
×
3
a
3
2
c
m
1
M
3
a
3
2
c
m
2
M
e
,
F
t
=
0
0
0
0
0
0
0
0
0
F
13
t
F
23
t
T
,where(14)m1=diag⁡m1m1m1,m2=diag⁡m2m2m2,m3=diag⁡m3m3m3,a1=a2a3a1,a2=-a3a5-a4,a3=-a3-a2a1,a4=a3-a5a4,k1ζ=diag⁡k1xk1yk1z,k2ζ=diag⁡k2xk2yk2z,k3ζ=diag⁡k3xk3yk3z,c1ζ=diag⁡c1xc1yc1z,c2ζ=diag⁡c2xc2yc2z,c3ζ=diag⁡c3xc3yc3z,a1ζ=-a2-a3-a1T,a2ζ=a3a2-a1T,a3ζ=a3-a5a4T,a3′ζ=-a3a5a4T,F13t=a3r1T0J1+a3r3T1J3-e¨1t,F23t=a3r2T0J2+a3r3T1J3-e¨2t.M
i and M
e are the equivalent masses introduced to simplify the computational process, defined by(15)Mi=Jiri2,Me=J2J3J2r32+J3r22.Mesh dampingc
m
i associated with the mesh stiffness can be obtained by(16)cmi=2ξikmi1/Mi+1/M30.5,where ξ
i is the damping ratio (ξ
i
=
0.01).Nonlinear backlash function is described by a piecewise formulationf
(
t
,
b
n
), given by(17)ft,bn=δijt-bn,δijt>bn,0,-bn≤δijt≤bn,δijt+bn,δijt<-bn.To enhance the ability of convergence in the differential equation computation, we introduce the termsb
n
′ and ω
n for the nondimensionalization operation, where b
n
′
=
100
×
1
0
-
6 m represents a characteristic length and ω
n is the characteristic frequency used to define dimensionless time τ as well as other quantities. Here ω
n is defined by(18)ωn=kmj1M1+1M30.5.With these terms, (12) yields the dimensionless equations of motion and is rewritten in the following forms: (19)m1ωn2bn′∂2ζ^1∂t2+c1ζ^1ωnbn′∂ζ^1∂t+k1ζ^1bn′ζ^1+a1ζ^1cm1ωnbn′δ˙13+km1tbn′f^δ13=0,m2ωn2bn′∂2ζ^2∂t2+c2ζ^2ωnbn′∂ζ^2∂t+k2ζ^2bn′ζ^2+a2ζ^2cm2ωnbn′δ˙23+km2tbn′f^δ23=0,m3ωn2bn′∂2ζ^3∂t2+c3ζ^3ωnbn′∂ζ^3∂t+k3ζ^3bn′ζ^3+a3ζ^3cm1ωnbn′δ˙13+km1tbn′f^δ13+a3′ζ^3cm2ωnbn′δ˙23+km2tbn′f^δ23=0,δ¨13=-a2X¨1-a3Y¨1-a1Z¨1+a3X¨3-a5Y¨3+a4Z¨3-2a32ξ1δ˙13-a32k^1ωn2Mef^δ13-2a32ξ2δ˙23-a32k^2ωn2Mef^δ23+a3r1T0ωn2J1bn′+a3r3T1ωn2J3bn′-e^¨1τωn,δ¨23=a3X¨2+a2Y¨2-a1Z¨2-a3X¨3+a5Y¨3+a4Z¨3-2a32ξ2δ˙23-a32k^2ωn2Mef^δ23-2a32ξ1δ˙13-a32k^1ωn2Mef^δ13+a3r2T0ωn2J2bn′+a3r3T1ωn2J3bn′-e^¨2τωn,where ζ
^
i
=
X
i
,
Y
i
,
Z
i
T, i
=
1,2
,
3.Other essential nondimensional parameters used are defined as follows:(20)τ=ωnt,ω=Ωωn,b=bnbn′,k^m1τ=km1t,k^m2τ=km2t,c^m1τ=cm1t,c^m2τ=cm2t,ξ1=c^m12ωnMe,ξ2=c^m22ωnMe,ζ^τ=ζtbn′,ζ˙^τ=ζ˙tωnbn′,ζ¨^τ=ζ¨tωn2bn′,f^τ,bn′=ft,bnbn′,e^¨iτ=e¨itωn2bn′,M^=M,K^=Kωn2,C^=Cωn,F^=Fωn2bn′.
^  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,λ
˙
1
=
λ
2, λ
˙
3
=
λ
4, λ
˙
5
=
λ
6, Δ
˙
1
=
Δ
2, and Δ
˙
3
=
Δ
4
, to make the system equations become a first-order differential equation set, given by(21)λ˙2=f1λ1,λ2,Δ1,Δ2λ˙4=f2λ3,λ4,Δ3,Δ4λ˙6=f3λ5,λ6,Δ1,Δ2,Δ3,Δ4Δ˙2=f4λ1,λ2,λ3,λ4,λ5,λ6,Δ1,Δ2,Δ3,Δ4Δ˙4=f5λ1,λ2,λ3,λ4,λ5,λ6,Δ1,Δ2,Δ3,Δ4,where λ
i
=
x
i
,
y
i
,
z
i, i
=
(
1,2
,
3
,
…
,
6
), represent the translational displacement vector and Δ
i (i
=
1,2
,
3,4) represent the relative displacements.


## 3. Bifurcation Behaviors

For (12), as a nonautonomous equation set with multiple parameters, it can be generally described as follows:(22)λ˙=Jλ,μ,τλ,λτ0=λ0,λ∈RN,μ∈RK,τ∈R,where J
(
λ
,
μ
,
τ
) is a Jacobi matrix, μ stands for the excitation source in the system, and λ
0 represents an initial condition.


### 3.1. Lyapunov Exponent

Introduce two close initial valuesλ
5
δ
13
(
0
)
,
δ
˙
13
(
0
) and λ
5
′
δ
13
′
(
0
)
,
δ
˙
13
′
(
0
), and define them by L
0
=
λ
5
-
λ
5
′
=
0.0001, where L
0 stands for the distance of those two state vectors. After several steps of computation, the initial value would become λ
5
(
δ
13
(
t
)
,
δ
˙
13
(
t
)
) and λ
5
′
δ
13
′
(
t
)
,
δ
˙
13
′
(
t
), and the separation of those two orbits can be evaluated through [24, 25](23)Lt=Lλ5,λ5′=sqrt⁡∑λ5-λ5′2.While (23) reaches the steady state, the state variable would develop into λ
5
(
n
τ
) and λ
5
′
(
n
τ
) at the time step of (
n
)th, where τ is the time step in Gill’s computation, n is the current number of iteration, and the distance at this moment is(24)Ln=lnτ=λ5nτ-λ5′nτ.With the iteration of the state points, the result on the orbits at the next step(
n
+
1
)th would be(25)λ5n+1τ=fτλ5nτ,λ5′n+1τ=fτλ5′nτ,where f
τ represent the mapping formula under parametric excitation. Combining (24) and (25) the trajectory separation at (
n
+
1
)th step will be(26)Ln+1=fτλ5nτ-fτλ5′nτ.One can find a start pointλ
5
′ that will keep a constant distance L
0 away from λ
5
[
(
n
+
1
)
τ
] which avoids the divergence of L
n
+
1 during numerical procedure. After a certain number of steps, a series of Lyapunov exponents are obtained as follows:(27)1τln⁡L1L0,1τln⁡L2L0,…,1τln⁡LiL0,where i
=
1,2
,
…
,
∞.The largest Lyapunov exponentλ
1 in terms of (27) can be obtained by(28)λ1=limn→∞⁡1nτ∑i=1nln⁡LiL0.


### 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 μ
j
=
u
j
T, where u
j is one of the concerned parameters, in the global domain with the forms of(29)μj=⋃j=1Nujuj=∏j=1nμj,μj+1,where μ
j is a unit cell in domain space.The centre of the cell is used in computing. For anm-dimensional domain, the cell identification can be discretized as(30)
Π
j
i
=
1
2
μ
1
,
j
+
μ
1
,
j
+
1
,
1
2
μ
2
,
j
+
μ
2
,
j
+
1
,
…
,
1
2
μ
i
,
j
+
μ
i
,
j
+
1
T
i
=
1,2
,
…
,
m
,where i refers to the quantities of interest and j refers to the jth 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 thejth cell is transformed as(31)Πj2=12μ1,j+μ1,j+1,12μ2,j+μ2,j+1T.Under the periodical excitation, the effect with regard toμ
j has the following relation to the system response:(32)Xn+1=Xn+∫t0t0+Tfμj,tdt.Data is only captured after the computation runs into steady state, according to (32), and supposing the system is in thePeriod-1 state, it can be described by X
i
+
1
=
C
(
X
i
,
μ
j
), if it is a multiple periodic motion, and the mapping route of term μ
j would be (33)Xi,μj⟶CXi,μj⟶C2Xi,μj⟶⋯⟶CkXi,μjk∈R.Ifk
→
∞, (33) is still not satisfied, which means that the cell of μ
j is in a chaotic state.In a dynamic regionΣ which contains n
×
n cells, the whole domain may consist of the periodic, quasi-periodic, or chaotic state due to the analytic solution. If the cells μ
j demonstrate being k-periodic motion, such cells can be described by (34)∑k=⋃j=1nμjk=1,2,…,m.The union of quasi-periodic cell is(35)∑q=⋃j=1nμjq=1,2,…,k.Similarly, it yields the chaotic domain(36)∑c=⋃j=1nμjifk>nT.The whole domain plane can be assembled by(37)∑=⋃k=1m∑k+⋃q=1k∑q+∑c.


### 3.3. Bifurcation Evolution

With respect to the bifurcation investigation, the key parameteru
j is set in a certain range and varies with time τ; accordingly a series of state vectors can be obtained via X
j
(
u
j
,
τ
), j
=
1,2
,
…
,
n. Here, the bifurcation diagram is considered using deflection displacement δ
13 versus variable u
j which can be plotted by (
u
j
,
δ
13
). Figure 2 demonstrates the threshold of bifurcation as well as the development process from a simple motion to chaos with the changes of E
n. At first E
n
=
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 thePeriod-2 orbit be split as two sections before coming intoPeriod-1 motion; meanwhile, the instability region emerges from ω
>
0.954. When E
n
=
28.7 and E
n
=
29, thePeriod-2 motion over the range ω
∈
[
0.935
,
0.951
] becomes more and more degenerate and accompanied by further and further strengthening of the unstable vibration at ω
>
0.951. For E
n
=
30 and E
n
=
32, the right side of thePeriod-2 orbit has disappeared, giving way to stablePeriod-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 E
n
=
42, we see the initial right branches ofPeriod-2 which have been totally replaced byPeriod-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 ofPeriod-9,Period-11, andPeriod-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.Figure 2
Bifurcation evolution procedure in terms of the mesh frequencyω
∈
[
0.9
,
1.08
].
(a)
E
n
=
28
(b)
E
n
=
28.7
(c)
E
n
=
29
(d)
E
n
=
30
(e)
E
n
=
32
(f)
E
n
=
42Displacement bifurcation versus dimensionless frequency has been shown in Figure3(a) based on ω
∈
[
1.68
,
2.21
]. One can find two significant skips that appear at ω
=
1.927 and ω
=
2.112. We note that the first one also breaks the route ofPeriod-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.Figure 3
Bifurcation behavior with respect to the mesh frequencyω
∈
[
1.68
,
2.21
].
(a)
Bifurcation diagram
(b)
Largest Lyapunov exponent


### 3.4. Domain Structure Analysis

Discretize the parameter space to be a two-dimensional plane withc
n
⊗
c
n, where c
n denotes the nth cell. The mapping initial condition is set as λ
5
=
0, letting ε
1
=
1
e
-
2 and ε
2
=
1
e
-
4 to be the computing precision to distinguish the target cell unit and identify the dynamic periodicity; for instance, if λ
5
(
i
)
-
λ
5
(
i
+
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 Figure4(a) the whole plane is divided into 22500 rectangular cells with uniform size of intervals of 0.008 and 0.0008 for E
n and b
n, respectively. Clearly, the global attraction is constructed with three distinct domains of P
1, P
4, andChaos, where P represents the periodic domains; the numbers 1,4 refer to the motion characteristics. It is obvious that the window ofPeriod-4 orbit is filled inside the chaotic domain. With increasing the transmission error E
n, 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 E
n
⊗
b
n plane in the case b
n
=
0.1 and E
n
∈
[
0.3
,
1.6
], 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 P
1 show that backlash plays a small role on the bifurcation adjustment when E
n
<
0.4. Another parameter domain plane E
n
⊗
m
k with [
0.3
,
1.6
]
⊗
[
0.05
,
0.4
] was conducted in Figure 4(b) and received the similar attractions with P
1, P
4, andChaos. Bifurcation diagram in Figure 4(d) indicates a good agreement with the parameter domain analysis.Figure 4
Domain structure analysis: (a), (b) parameter planes; (c), (d) bifurcation diagrams.
(a)
E
n
⊗
b
n
∈
[
0.3
,
1.6
]
⊗
[
0.04
,
0.18
]
(b)
E
n
⊗
m
k
∈
[
0.3
,
1.6
]
⊗
[
0.05
,
0.4
]
(c)
E
n
∈
[
0.3
,
1.6
], b
n
=
0.1
(d)
m
k
∈
[
0.05
,
0.4
], E
n
=
1.2The results in Figure5 demonstrate the complicated bifurcation evolutions in the bounded domain. The dynamic group mainly consists of P
1, P
2, and P
3 as well asChaos, 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 E
n
⊗
ω plane,Period-2 cells have a spatial domination over the system response. During the changes of decreasing E
n, thePeriod-1 cells disappeared firstly and thePeriod-3 orbit moves upwards along the C edge on the left side ofPeriod-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 fromPeriod-1 motion to an irregular region and then to chaotic domain, and chaos appears followingPeriod-3 cells.Figure 5
Parametric plane for the case.E
n
⊗
ω
∈
[
1.2
,
2.5
]
⊗
[
0.74
,
0.825
]. T
0
=
400, ξ
=
0.01, b
=
0.093, and m
k
=
0.08.In Figure6(a), the occurrence of bifurcation motion is synchronous with largest Lyapunov exponent in the case ω
=
0.788, and the system vibration goes through stablePeriod-1 andPeriod-2 orbit to the chaotic vibration. Bifurcation diagram and the largest LE both discover that a sudden jumping occurs at E
n
=
1.45 just happening atPeriod-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.Figure 6
Verification of the domain structure with respect to transmission error.
(a)
Bifurcation diagram withω
=
0.788
(b)
Largest Lyapunov exponentFor Figure7, the CMM solutions are depicted employing the damping ratio of ξ
∈
[
0.05
,
0.41
] and, respectively, versus parameters of backlash b
n, mesh error E
n, stiffness factor m
k, and input torque T
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 P
2, P
4, P
8, and P
16; besides, some P
32 points are also evident at the edge before chaos emerges. In Figures 7(a)–7(c) afterPeriod-4 mapping region, small quantity ofPeriod-8 andPeriod-16 cells can be seen. In ξ
⊗
T
1 plane, very fewPeriod-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.Figure 7
Domain structure analysis in two-dimensional parameter planes.
(a)
ξ
⊗
b
∈
[
0.05
,
0.45
]
⊗
[
0.5
,
0.8
]
(b)
ξ
⊗
E
n
∈
[
0.05
,
0.45
]
⊗
[
0.25
,
0.4
]
(c)
ξ
⊗
m
k
∈
[
0.05
,
0.45
]
⊗
[
0.05
,
0.5
]
(d)
ξ
⊗
T
1
∈
[
0.05
,
0.45
]
⊗
[
100
,
230
]The condition in Figure8 is pictured with 300 steady periods for a fixed ξ. With the damping ratio decreasing, it mainly illustrates the changes starting from aPeriod-2
m response (m: integer) to a doubled periodic (m
+
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.Figure 8
Bifurcation evolution with respect to the damping ratio.
(a)
Bifurcation diagram withb
n
=
0.73
(b)
Largest Lyapunov exponent


## 3.1. Lyapunov Exponent

Introduce two close initial valuesλ
5
δ
13
(
0
)
,
δ
˙
13
(
0
) and λ
5
′
δ
13
′
(
0
)
,
δ
˙
13
′
(
0
), and define them by L
0
=
λ
5
-
λ
5
′
=
0.0001, where L
0 stands for the distance of those two state vectors. After several steps of computation, the initial value would become λ
5
(
δ
13
(
t
)
,
δ
˙
13
(
t
)
) and λ
5
′
δ
13
′
(
t
)
,
δ
˙
13
′
(
t
), and the separation of those two orbits can be evaluated through [24, 25](23)Lt=Lλ5,λ5′=sqrt⁡∑λ5-λ5′2.While (23) reaches the steady state, the state variable would develop into λ
5
(
n
τ
) and λ
5
′
(
n
τ
) at the time step of (
n
)th, where τ is the time step in Gill’s computation, n is the current number of iteration, and the distance at this moment is(24)Ln=lnτ=λ5nτ-λ5′nτ.With the iteration of the state points, the result on the orbits at the next step(
n
+
1
)th would be(25)λ5n+1τ=fτλ5nτ,λ5′n+1τ=fτλ5′nτ,where f
τ represent the mapping formula under parametric excitation. Combining (24) and (25) the trajectory separation at (
n
+
1
)th step will be(26)Ln+1=fτλ5nτ-fτλ5′nτ.One can find a start pointλ
5
′ that will keep a constant distance L
0 away from λ
5
[
(
n
+
1
)
τ
] which avoids the divergence of L
n
+
1 during numerical procedure. After a certain number of steps, a series of Lyapunov exponents are obtained as follows:(27)1τln⁡L1L0,1τln⁡L2L0,…,1τln⁡LiL0,where i
=
1,2
,
…
,
∞.The largest Lyapunov exponentλ
1 in terms of (27) can be obtained by(28)λ1=limn→∞⁡1nτ∑i=1nln⁡LiL0.


## 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 μ
j
=
u
j
T, where u
j is one of the concerned parameters, in the global domain with the forms of(29)μj=⋃j=1Nujuj=∏j=1nμj,μj+1,where μ
j is a unit cell in domain space.The centre of the cell is used in computing. For anm-dimensional domain, the cell identification can be discretized as(30)
Π
j
i
=
1
2
μ
1
,
j
+
μ
1
,
j
+
1
,
1
2
μ
2
,
j
+
μ
2
,
j
+
1
,
…
,
1
2
μ
i
,
j
+
μ
i
,
j
+
1
T
i
=
1,2
,
…
,
m
,where i refers to the quantities of interest and j refers to the jth 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 thejth cell is transformed as(31)Πj2=12μ1,j+μ1,j+1,12μ2,j+μ2,j+1T.Under the periodical excitation, the effect with regard toμ
j has the following relation to the system response:(32)Xn+1=Xn+∫t0t0+Tfμj,tdt.Data is only captured after the computation runs into steady state, according to (32), and supposing the system is in thePeriod-1 state, it can be described by X
i
+
1
=
C
(
X
i
,
μ
j
), if it is a multiple periodic motion, and the mapping route of term μ
j would be (33)Xi,μj⟶CXi,μj⟶C2Xi,μj⟶⋯⟶CkXi,μjk∈R.Ifk
→
∞, (33) is still not satisfied, which means that the cell of μ
j is in a chaotic state.In a dynamic regionΣ which contains n
×
n cells, the whole domain may consist of the periodic, quasi-periodic, or chaotic state due to the analytic solution. If the cells μ
j demonstrate being k-periodic motion, such cells can be described by (34)∑k=⋃j=1nμjk=1,2,…,m.The union of quasi-periodic cell is(35)∑q=⋃j=1nμjq=1,2,…,k.Similarly, it yields the chaotic domain(36)∑c=⋃j=1nμjifk>nT.The whole domain plane can be assembled by(37)∑=⋃k=1m∑k+⋃q=1k∑q+∑c.


## 3.3. Bifurcation Evolution

With respect to the bifurcation investigation, the key parameteru
j is set in a certain range and varies with time τ; accordingly a series of state vectors can be obtained via X
j
(
u
j
,
τ
), j
=
1,2
,
…
,
n. Here, the bifurcation diagram is considered using deflection displacement δ
13 versus variable u
j which can be plotted by (
u
j
,
δ
13
). Figure 2 demonstrates the threshold of bifurcation as well as the development process from a simple motion to chaos with the changes of E
n. At first E
n
=
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 thePeriod-2 orbit be split as two sections before coming intoPeriod-1 motion; meanwhile, the instability region emerges from ω
>
0.954. When E
n
=
28.7 and E
n
=
29, thePeriod-2 motion over the range ω
∈
[
0.935
,
0.951
] becomes more and more degenerate and accompanied by further and further strengthening of the unstable vibration at ω
>
0.951. For E
n
=
30 and E
n
=
32, the right side of thePeriod-2 orbit has disappeared, giving way to stablePeriod-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 E
n
=
42, we see the initial right branches ofPeriod-2 which have been totally replaced byPeriod-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 ofPeriod-9,Period-11, andPeriod-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.Figure 2
Bifurcation evolution procedure in terms of the mesh frequencyω
∈
[
0.9
,
1.08
].
(a)
E
n
=
28
(b)
E
n
=
28.7
(c)
E
n
=
29
(d)
E
n
=
30
(e)
E
n
=
32
(f)
E
n
=
42Displacement bifurcation versus dimensionless frequency has been shown in Figure3(a) based on ω
∈
[
1.68
,
2.21
]. One can find two significant skips that appear at ω
=
1.927 and ω
=
2.112. We note that the first one also breaks the route ofPeriod-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.Figure 3
Bifurcation behavior with respect to the mesh frequencyω
∈
[
1.68
,
2.21
].
(a)
Bifurcation diagram
(b)
Largest Lyapunov exponent


## 3.4. Domain Structure Analysis

Discretize the parameter space to be a two-dimensional plane withc
n
⊗
c
n, where c
n denotes the nth cell. The mapping initial condition is set as λ
5
=
0, letting ε
1
=
1
e
-
2 and ε
2
=
1
e
-
4 to be the computing precision to distinguish the target cell unit and identify the dynamic periodicity; for instance, if λ
5
(
i
)
-
λ
5
(
i
+
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 Figure4(a) the whole plane is divided into 22500 rectangular cells with uniform size of intervals of 0.008 and 0.0008 for E
n and b
n, respectively. Clearly, the global attraction is constructed with three distinct domains of P
1, P
4, andChaos, where P represents the periodic domains; the numbers 1,4 refer to the motion characteristics. It is obvious that the window ofPeriod-4 orbit is filled inside the chaotic domain. With increasing the transmission error E
n, 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 E
n
⊗
b
n plane in the case b
n
=
0.1 and E
n
∈
[
0.3
,
1.6
], 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 P
1 show that backlash plays a small role on the bifurcation adjustment when E
n
<
0.4. Another parameter domain plane E
n
⊗
m
k with [
0.3
,
1.6
]
⊗
[
0.05
,
0.4
] was conducted in Figure 4(b) and received the similar attractions with P
1, P
4, andChaos. Bifurcation diagram in Figure 4(d) indicates a good agreement with the parameter domain analysis.Figure 4
Domain structure analysis: (a), (b) parameter planes; (c), (d) bifurcation diagrams.
(a)
E
n
⊗
b
n
∈
[
0.3
,
1.6
]
⊗
[
0.04
,
0.18
]
(b)
E
n
⊗
m
k
∈
[
0.3
,
1.6
]
⊗
[
0.05
,
0.4
]
(c)
E
n
∈
[
0.3
,
1.6
], b
n
=
0.1
(d)
m
k
∈
[
0.05
,
0.4
], E
n
=
1.2The results in Figure5 demonstrate the complicated bifurcation evolutions in the bounded domain. The dynamic group mainly consists of P
1, P
2, and P
3 as well asChaos, 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 E
n
⊗
ω plane,Period-2 cells have a spatial domination over the system response. During the changes of decreasing E
n, thePeriod-1 cells disappeared firstly and thePeriod-3 orbit moves upwards along the C edge on the left side ofPeriod-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 fromPeriod-1 motion to an irregular region and then to chaotic domain, and chaos appears followingPeriod-3 cells.Figure 5
Parametric plane for the case.E
n
⊗
ω
∈
[
1.2
,
2.5
]
⊗
[
0.74
,
0.825
]. T
0
=
400, ξ
=
0.01, b
=
0.093, and m
k
=
0.08.In Figure6(a), the occurrence of bifurcation motion is synchronous with largest Lyapunov exponent in the case ω
=
0.788, and the system vibration goes through stablePeriod-1 andPeriod-2 orbit to the chaotic vibration. Bifurcation diagram and the largest LE both discover that a sudden jumping occurs at E
n
=
1.45 just happening atPeriod-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.Figure 6
Verification of the domain structure with respect to transmission error.
(a)
Bifurcation diagram withω
=
0.788
(b)
Largest Lyapunov exponentFor Figure7, the CMM solutions are depicted employing the damping ratio of ξ
∈
[
0.05
,
0.41
] and, respectively, versus parameters of backlash b
n, mesh error E
n, stiffness factor m
k, and input torque T
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 P
2, P
4, P
8, and P
16; besides, some P
32 points are also evident at the edge before chaos emerges. In Figures 7(a)–7(c) afterPeriod-4 mapping region, small quantity ofPeriod-8 andPeriod-16 cells can be seen. In ξ
⊗
T
1 plane, very fewPeriod-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.Figure 7
Domain structure analysis in two-dimensional parameter planes.
(a)
ξ
⊗
b
∈
[
0.05
,
0.45
]
⊗
[
0.5
,
0.8
]
(b)
ξ
⊗
E
n
∈
[
0.05
,
0.45
]
⊗
[
0.25
,
0.4
]
(c)
ξ
⊗
m
k
∈
[
0.05
,
0.45
]
⊗
[
0.05
,
0.5
]
(d)
ξ
⊗
T
1
∈
[
0.05
,
0.45
]
⊗
[
100
,
230
]The condition in Figure8 is pictured with 300 steady periods for a fixed ξ. With the damping ratio decreasing, it mainly illustrates the changes starting from aPeriod-2
m response (m: integer) to a doubled periodic (m
+
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.Figure 8
Bifurcation evolution with respect to the damping ratio.
(a)
Bifurcation diagram withb
n
=
0.73
(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 λ
1
∈
[
-
0.021
,
0.038
]; when ω
=
1.93, the discontinuous jump breaks thePeriod-2 orbit and at the point of ω
=
2.12, the jump switchesPeriod-2 state into chaos suddenly. As the transmission error changes from 28 to 42, bifurcation evolution was exhibited.
(2)
InE
n
⊗
b
n and E
n
⊗
m
k parameter spaces, three dynamic regions includingPeriod-1,Period-4, andChaos were exhibited, and gear dynamic responses are sensitive to the boundary parameters of the cell groups.
(3)
Inξ
⊗
b
n, ξ
⊗
E
n, ξ
⊗
m
k, and ξ
⊗
T
1 planes, the parameter domain structure reveals the global dynamic behavior distributions, with periodic regions of P
e
r
i
o
d
-
2
,
P
e
r
i
o
d
-
4
,
P
e
r
i
o
d
-
8
,
…,  andPeriod-2
m (m: integer), leading to chaos, and large damping ratio contributes to a steady dynamic response.



---

*Source: 3738508-2016-05-10.xml*