# Bifurcation and Chaos Prediction in Nonlinear Gear Systems

Authors: Anooshirvan Farshidianfar; Amin Saghafi
Journal: Shock and Vibration
 (2014)

---

## Abstract

The homoclinic bifurcation and transition to chaos in gear systems are studied both analytically and numerically. Applying Melnikov analytical method, the threshold values for the occurrence of chaotic motion are obtained. The influence of system parameters on the character of vibration is studied. The numerical simulation of the system including bifurcation diagram, phase plane portraits, Fourier spectra, and time histories is considered to confirm the analytical predictions for the occurrence of homoclinic bifurcation and chaos in nonlinear gear systems.


---

## Body

## 1. Introduction

Gear systems are known as one of the important sources of noise and vibration in industrial rotating machinery and power transmission systems. With the development of nonlinear dynamics, the nonlinear characteristics in gear systems, such as stability, periodic responses, bifurcations, and chaotic behaviors, have become the most interesting research areas. For instance, Sato et al. [1] established a nonlinear model of gear system with the time dependence of tooth stiffness and backlash. They investigated the bifurcation and chaotic phenomena by using a shooting method. Padmanabhan and Singh [2] analyzed a 2-d.o.f system with clearance nonlinearity. They used a numerical simulation to demonstrate the existence of periodic and chaotic behaviors. Kahraman and Blankenship [3] performed experiments on a spur gear pair and observed various nonlinear phenomena including period doubling and chaos. Theodossiades and Natsiavas [4] investigated the motor-driven gear-pair systems and found periodic and chaotic behaviors in this system. Also, Chang-Jian and Chang [5] investigated the dynamic responses of a single-degree-of freedom spur gear system with and without nonlinear suspension and found the bifurcation and chaotic dynamics in this system.From the above-mentioned references, one finds that chaotic phenomena have been found in nonlinear gear systems. Analytical [2], numerical [1, 4, 5], and experimental methods [3] were used by different authors to study chaos for spur gear systems. Due to the complexity of the gear systems and also the difficulty and limitation of the analytical methods, the numerical method was commonly used to analyze the nonlinear dynamics in gear systems, but they cannot provide any analytical expression of the solutions. Although significant amount of research has been devoted to nonlinear gear dynamic, few attempts have been made to investigate analytical solution of the system. Melnikov analysis is one of the few analytical methods to provide an approximate criterion for the occurrence of hetero/homoclinic bifurcation and chaos in nonlinear systems. According to this theory, the existence of transversal intersection of stable and unstable manifolds of saddle fixed points implies the occurrence of the chaos [6–8]. Thus, the current paper focuses on the study of the homoclinic bifurcation and chaos for the gear system by means of Melnikov analysis.The organization of the paper is as follows. In Section2, a model for the vibration of gear transmission system is presented. In Section 3, the conditions for chaotic behaviour in terms of homoclinic bifurcation are given by using Melnikov analysis. Finally, in Section 4, the threshold values for the occurrence of chaotic motion are obtained. Some numerical simulation of the system including bifurcation diagrams, plane phase portraits, Fourier spectra, and time histories is used to confirm the analytical predictions and show the transition to chaotic motion.


## 2. Problem Formulation

The mechanical model of a spur gear pair system investigated in the present study is shown in Figure1(a). In this model, the gear mesh is modeled as a pair of rigid disks connected by a spring damper set along the line of action. The backlash function is usually used to represent gear clearances. Under these assumptions, the equations of motion can be expressed as(1a)
I
a
d
2
θ
a
d
t
-
2
+
c
(
r
a
d
θ
a
d
t
-
-
r
b
d
θ
b
d
t
-
)
r
a
+
r
a
K
h
(
r
a
θ
a
-
r
b
θ
b
)
=
T
a
,
(1b)
I
b
d
2
θ
b
d
t
-
2
-
c
(
r
a
d
θ
a
d
t
-
-
r
b
d
θ
b
d
t
-
)
r
b
-
r
b
K
h
(
r
a
θ
a
-
r
b
θ
b
)
=
-
T
b
.In these equations, θ
a and θ
b are the torsional displacements of gears (a) and (b). I
a and I
b are the mass moments of inertia and r
a and r
b are the base circle radius of the gears. K and c represent the stiffness and linear viscous damping of the gear mesh. T
a and T
b are external torques acting on the driver and driven gear, respectively. An input torque T
a is applied to the driver gear (a) rotating at Ω
a and the mean braking torque T
b on the driven gear with angular velocity Ω
b. The excitation torque T
a fluctuates significantly between low and high values. Therefore, the T
a can be decomposed into average torque transmitted through the gear pair T
m
a and the fluctuating external torque excitation T
p
(
t
-
) parts. Such excitations are typically at low frequencies Ω which are the first few multiples of the input shaft frequency. Also, output torque fluctuation will be neglected; that is, T
b
(
t
-
)
=
T
m
b. Express T
a via Fourier series as [9]
(2)
T
a
(
t
-
)
=
T
m
a
+
∑
r
=
1
∞
T
p
r
cos
⁡
(
r
Ω
t
-
+
φ
r
)
.
So, (1a) and (1b) are simplified into (3) by introducing a new variable x
-
=
r
a
θ
a
-
r
b
θ
b:
(3)
m
x
-
¨
+
c
x
-
˙
+
k
h
(
x
-
)
=
F
-
m
+
∑
r
=
1
∞
F
-
p
r
cos
⁡
(
r
Ω
t
-
+
φ
r
)
with
(4)
h
(
x
-
)
=
{
x
-
-
(
1
-
α
)
d
d
<
x
-
α
x
-
-
d
≤
x
-
≤
d
,
x
-
+
(
1
-
α
)
d
d
<
-
x
-
m
=
I
a
I
b
I
b
r
a
2
+
I
a
r
b
2
,
F
-
p
r
=
m
(
r
a
I
a
)
T
p
r
,
F
-
m
=
m
(
T
m
a
r
a
I
a
+
T
m
b
r
b
I
b
)
,
Here m is the equivalent mass representing the total inertia of the gear pair, F
-
m is the average force transmitted through the gear pair, and F
-
p
(
t
-
) is the fluctuating force related to the input torque excitation. 2
d represent the backlash. h
(
x
-
) is the nonlinear displacement function due to backlash. It is a stepwise linear function. A nondimensional form of the above equation is obtained by letting
(5)
Ω
n
=
k
m
,
t
=
Ω
n
t
-
,
x
=
x
-
d
c
,
μ
^
=
c
m
Ω
n
,
F
^
m
=
F
-
m
d
c
k
,
F
^
p
r
=
F
-
p
r
d
c
k
,
ω
=
Ω
Ω
n
,
where d
c is characteristic length. By choosing the parameters α
=
0 and d
/
d
c
=
0.1745 [9], the function h
(
x
) is approximated with a third-order polynomial as h
(
x
)
=
-
0.1422
x
+
5.8093
x
3 (Figure 1(b)). Substituting h
(
x
), the nondimensional form of (3) could be rewritten as
(6)
x
¨
+
μ
^
x
˙
+
(
-
0.1422
x
+
5.8093
x
3
)
=
F
^
m
+
∑
r
=
1
∞
F
^
p
r
cos
⁡
(
ω
t
+
φ
r
)
.
This equation considers the effect of backlash and so it is the equation which would be studied in this paper.(a) A gear pair model and (b) the approximation ofh
(
x
) based on a third-order polynomial.
(a)
(b)


## 3. Global Bifurcation and Chaos in Gear Model

The study of homoclinic bifurcation that enables predicting the chaotic behaviors of nonlinear systems is well done by the Melnikov theory. The Melnikov method is one of the few analytical methods to study the global bifurcation of the system and gives a procedure for analyzing and estimating when a chaotic behavior of a nonlinear system is expected. According to this theory, the existence of transversal intersection of stable and unstable manifolds of homoclinic orbits implies the existence of the chaotic dynamic. In order to apply this technique and to carry out this study, the average force, the excitation term, and the damping are considered as small perturbations for the Hamiltonian system. Scalingμ
^
=
ε
μ, F
^
m
=
ε
f
m, and F
^
p
r
=
ε
f
p
r, (6) can be rewritten as
(7)
x
˙
=
y
,
y
˙
=
-
ε
μ
x
˙
+
(
0.1422
x
-
5.8093
x
3
)
+
ε
(
f
m
+
∑
r
=
1
∞
f
p
r
cos
⁡
⁡
(
ω
t
+
φ
r
)
)
,
where ε is a small parameter so that the damping, forcing, and excitation terms are perturbations to the unperturbed system.


### 3.1. Analysis of the Unperturbed Model

In this subsection, the homoclinic orbits of the unperturbed system are derived. For the unperturbed system, whenε
=
0, (7) is simplified to
(8)
x
˙
=
y
,
y
˙
=
(
b
x
-
c
x
3
)
,
where b
=
0.1422 and c
=
5.8093. The system of (8) is a Hamiltonian system with a potential function as
(9)
U
(
x
)
=
-
b
x
2
2
+
c
x
4
4
.
The Hamiltonian function of the unperturbed system is given by
(10)
H
(
x
,
y
)
=
y
2
2
+
U
(
x
)
=
1
2
y
2
-
b
2
x
2
+
c
4
x
4
.
The potential function and the Hamiltonian function of the unperturbed system are shown in Figures 2(a) and 2(b), respectively. The unperturbed system has three fixed points: (
0,0
) and (
±
b
/
c
,
0
). From the linear stability analysis, (
b
/
c
,
0
) and (
-
b
/
c
,
0
) are centers, and (
0,0
) is a saddle point. The saddle point is connected to itself by two homoclinic orbits. In order to obtain homoclinic function, transforming (10) for a constant energy, chosen as zero, the homoclinic trajectories are obtained as
(11)
(
x
h
±
(
τ
)
,
y
h
±
(
τ
)
)
=
(
±
2
b
c
sec
h
(
b
(
τ
)
)
,
(
x
h
±
(
τ
)
,
y
h
±
(
τ
)
)
=
∓
2
c
b
sec
h
(
b
(
τ
)
)
tanh
(
b
(
τ
)
)
)
,
where t
-
t
0
=
τ. Stable and unstable manifolds of homoclinic orbits (W
s
± and W
u
±) are shown in Figure 3.(a) Potential function and (b) Hamiltonian function of the unperturbed system.
(a)
(b)Figure 3
Phase portrait and homoclinic orbits of unperturbed system.


### 3.2. Melnikov Analysis for Gear Equation

In this subsection, the conditions for existence of the homoclinic bifurcation and chaos by using the Melnikov method are given. The Melnikov function measures the distance between the stable and unstable manifolds in the Poincare section. According to this theory, the existence of transversal intersection of stable and unstable manifolds of homoclinic orbits implies the existence of the chaotic dynamics. The generalized equation (see (7)) can be written in the vector form as
(12)
x
˙
=
p
1
(
x
,
y
)
+
ε
q
1
(
x
,
y
,
t
)
,
y
˙
=
p
2
(
x
,
y
)
+
ε
q
2
(
x
,
y
,
t
)
,
where
(13)
p
1
(
x
,
y
)
=
y
,
q
1
(
x
,
y
,
t
)
=
0
,
p
2
(
x
,
y
)
=
(
b
x
-
c
x
3
)
,
q
2
(
x
,
y
,
t
)
=
-
μ
x
˙
+
f
m
+
∑
r
=
1
∞
f
p
r
cos
⁡
(
ω
t
+
φ
r
)
.
Using the Melnikov function the distance between the stable and unstable manifolds is defined as [6]
(14)
M
(
t
0
)
=
∫
-
∞
+
∞
p
(
X
h
(
τ
)
)
∧
q
(
X
h
(
τ
)
,
t
)
×
exp
⁡
[
-
∫
0
T
trace
[
D
X
(
p
(
X
h
(
s
)
)
)
]
d
s
]
d
τ
,
where X
h
=
(
x
h
,
y
h
) represents homoclinic orbits, D
X is the partial derivative with respect to X
=
(
x
,
y
), and p
∧
q
=
p
1
q
2
-
p
2
q
1. Substituting (12) into (14), the Melnikov integral could be rewritten as
(15)
M
(
t
0
)
±
=
∫
-
∞
+
∞
[
∓
2
c
b
sec
h
(
b
τ
)
tanh
⁡
(
b
τ
)
]
⋯
[
-
μ
(
∓
2
c
b
sec
h
(
b
τ
)
tanh
(
b
τ
)
)
+
f
m
+
∑
r
=
1
∞
f
p
r
cos
⁡
⁡
(
ω
(
τ
+
t
0
)
+
φ
r
)
]
d
τ
.
Considering only the first harmonics term (r
=
1) and evaluating the integral, the Melnikov function can be written as
(16)
M
(
t
0
)
±
=
-
4
3
μ
(
b
)
2
c
b
±
2
c
f
p
π
ω
sec
h
(
π
ω
2
b
)
sin
⁡
(
ω
t
0
+
φ
r
)
=
A
±
D
sin
(
ω
t
0
+
φ
r
)
.
If the Melnikov function has a simple zero, then the stable and unstable manifolds intersect transversally, homoclinic bifurcation occurs and hence chaotic behavior is expected. Thus, the conditions for the occurrence of chaos are obtained as
(17)
|
A
|
≤
|
D
|
⟹
|
-
4
3
μ
(
b
)
2
c
b
|
≤
|
2
c
f
p
π
ω
sec
h
(
π
ω
2
b
)
|
.
This condition provides a domain on the parameter spaces where the system has transverse homoclinic orbits resulting in possible chaotic behavior.


## 3.1. Analysis of the Unperturbed Model

In this subsection, the homoclinic orbits of the unperturbed system are derived. For the unperturbed system, whenε
=
0, (7) is simplified to
(8)
x
˙
=
y
,
y
˙
=
(
b
x
-
c
x
3
)
,
where b
=
0.1422 and c
=
5.8093. The system of (8) is a Hamiltonian system with a potential function as
(9)
U
(
x
)
=
-
b
x
2
2
+
c
x
4
4
.
The Hamiltonian function of the unperturbed system is given by
(10)
H
(
x
,
y
)
=
y
2
2
+
U
(
x
)
=
1
2
y
2
-
b
2
x
2
+
c
4
x
4
.
The potential function and the Hamiltonian function of the unperturbed system are shown in Figures 2(a) and 2(b), respectively. The unperturbed system has three fixed points: (
0,0
) and (
±
b
/
c
,
0
). From the linear stability analysis, (
b
/
c
,
0
) and (
-
b
/
c
,
0
) are centers, and (
0,0
) is a saddle point. The saddle point is connected to itself by two homoclinic orbits. In order to obtain homoclinic function, transforming (10) for a constant energy, chosen as zero, the homoclinic trajectories are obtained as
(11)
(
x
h
±
(
τ
)
,
y
h
±
(
τ
)
)
=
(
±
2
b
c
sec
h
(
b
(
τ
)
)
,
(
x
h
±
(
τ
)
,
y
h
±
(
τ
)
)
=
∓
2
c
b
sec
h
(
b
(
τ
)
)
tanh
(
b
(
τ
)
)
)
,
where t
-
t
0
=
τ. Stable and unstable manifolds of homoclinic orbits (W
s
± and W
u
±) are shown in Figure 3.(a) Potential function and (b) Hamiltonian function of the unperturbed system.
(a)
(b)Figure 3
Phase portrait and homoclinic orbits of unperturbed system.


## 3.2. Melnikov Analysis for Gear Equation

In this subsection, the conditions for existence of the homoclinic bifurcation and chaos by using the Melnikov method are given. The Melnikov function measures the distance between the stable and unstable manifolds in the Poincare section. According to this theory, the existence of transversal intersection of stable and unstable manifolds of homoclinic orbits implies the existence of the chaotic dynamics. The generalized equation (see (7)) can be written in the vector form as
(12)
x
˙
=
p
1
(
x
,
y
)
+
ε
q
1
(
x
,
y
,
t
)
,
y
˙
=
p
2
(
x
,
y
)
+
ε
q
2
(
x
,
y
,
t
)
,
where
(13)
p
1
(
x
,
y
)
=
y
,
q
1
(
x
,
y
,
t
)
=
0
,
p
2
(
x
,
y
)
=
(
b
x
-
c
x
3
)
,
q
2
(
x
,
y
,
t
)
=
-
μ
x
˙
+
f
m
+
∑
r
=
1
∞
f
p
r
cos
⁡
(
ω
t
+
φ
r
)
.
Using the Melnikov function the distance between the stable and unstable manifolds is defined as [6]
(14)
M
(
t
0
)
=
∫
-
∞
+
∞
p
(
X
h
(
τ
)
)
∧
q
(
X
h
(
τ
)
,
t
)
×
exp
⁡
[
-
∫
0
T
trace
[
D
X
(
p
(
X
h
(
s
)
)
)
]
d
s
]
d
τ
,
where X
h
=
(
x
h
,
y
h
) represents homoclinic orbits, D
X is the partial derivative with respect to X
=
(
x
,
y
), and p
∧
q
=
p
1
q
2
-
p
2
q
1. Substituting (12) into (14), the Melnikov integral could be rewritten as
(15)
M
(
t
0
)
±
=
∫
-
∞
+
∞
[
∓
2
c
b
sec
h
(
b
τ
)
tanh
⁡
(
b
τ
)
]
⋯
[
-
μ
(
∓
2
c
b
sec
h
(
b
τ
)
tanh
(
b
τ
)
)
+
f
m
+
∑
r
=
1
∞
f
p
r
cos
⁡
⁡
(
ω
(
τ
+
t
0
)
+
φ
r
)
]
d
τ
.
Considering only the first harmonics term (r
=
1) and evaluating the integral, the Melnikov function can be written as
(16)
M
(
t
0
)
±
=
-
4
3
μ
(
b
)
2
c
b
±
2
c
f
p
π
ω
sec
h
(
π
ω
2
b
)
sin
⁡
(
ω
t
0
+
φ
r
)
=
A
±
D
sin
(
ω
t
0
+
φ
r
)
.
If the Melnikov function has a simple zero, then the stable and unstable manifolds intersect transversally, homoclinic bifurcation occurs and hence chaotic behavior is expected. Thus, the conditions for the occurrence of chaos are obtained as
(17)
|
A
|
≤
|
D
|
⟹
|
-
4
3
μ
(
b
)
2
c
b
|
≤
|
2
c
f
p
π
ω
sec
h
(
π
ω
2
b
)
|
.
This condition provides a domain on the parameter spaces where the system has transverse homoclinic orbits resulting in possible chaotic behavior.


## 4. Numerical Simulations

In this section, we give the numerical simulation to demonstrate the theoretical results from Melnikov analysis obtained in the previous section. The influence of control parameters on the character of vibration and frequency response is studied. Using (17) and choosing f
p as control parameter, the condition for transverse intersection of stable and unstable manifolds is obtained as
(18)
f
p
≥
2
2
3
μ
(
b
)
3
/
2
π
ω
c
cos
⁡
h
(
π
ω
2
b
)
.
This condition provides a domain on the parameter spaces where the system has transverse homoclinic orbits resulting in possible chaotic dynamics. Figure 4 depicts the threshold curves f
p for homoclinic bifurcation in the (
f
p
-
ω
) plane for μ
=
2. In the parameter region above the threshold curve the transverse intersections of stable and unstable manifolds occur. As it can be observed from this figure, the threshold value is f
p
=
0.43, for ω
=
1. Since Melnikov function can change its sign for f
p
>
0.43, chaos may occur.Figure 4
Homoclinic bifurcation curve in the (f
p
-
ω) plane.To confirm the analytical prediction, the bifurcation diagram has been depicted. The nonlinear equation (6) is integrated numerically using the fourth-order Runge-Kutta method. Figure 5 presents the bifurcation diagram of system using the f
p as a bifurcation parameter. The values of the parameters are fixed as ω
=
1, μ
=
2, f
m
=
0.05, and ε
=
0.1 and initial conditions as x
=
0.01 and x
˙
=
0.01. It can be observed that the gear system exhibits a 1T-periodic response at low values of f
p; that is, f
p
≤
0.35. However, as f
p is increased the period 1 solution bifurcates to a 2T-periodic solution through a period doubling bifurcation. With the increase of the parameter f
p, period doubling occurs, and a bifurcation cascade leading to chaos, which is consistent with the prediction of analytical chaos, was presented in Figure 4.Figure 5
Bifurcation diagram usingF
^
p (F
^
p
=
ε
f
p) as control parameter.For a better clarity, the transition to chaos through the numerical simulation is shown for three values off
p
=
0.3, f
p
=
0.38, and f
p
=
0.45 (see Figures 4 and 5). Point 1 (f
p
=
0.3, ω
=
1) and point 2 (f
p
=
0.38, ω
=
1), lying out of chaos, and point 3 (f
p
=
0.45, ω
=
1), lying in the chaotic area, can be observed. Point 1 corresponds to the point situated below threshold values and is associated with the periodic motion. The numerical simulations of (6) are performed for the parameter sets associated with this point. The time history, phase plane diagram, and Fourier spectra are shown in Figure 6. These figures show the periodic response and confirm that the prediction of a periodic motion was obtained by Melnikov theory.(a) Time history, (b) Fourier spectra, and (c) phase plane, forf
p
=
0.3.
(a)
(b)
(c)In the next computational step,f
p
=
0.38 (point 2) is considered. It corresponds to the point situated below threshold values (Figure 4) and is associated with the 2T-periodic motion (Figure 5). The Melnikov function does not change its sign; thus, a periodic motion occurs. Figure 7 illustrates the time history, phase plane, and Fourier spectra for the parameter f
p
=
0.38. Observe that for this point the 2T-periodic motion occurs.(a) Time history, (b) Fourier spectra, and (c) phase plane, forf
p
=
0.38.
(a)
(b)
(c)Finally, the numerical analysis is carried out for the parameterf
p
=
0.45. The Melnikov function changes its sign for f
p
=
0.45 and corresponds to the chaotic motion. The chaotic trajectories are obtained by numerical integration. The time history, phase plane, and Fourier spectra corresponding to the chaotic response are presented in Figure 8. Chaotic behavior is clearly visible. The numerical computations confirm the prediction of analytical chaos for applied value of f
p.(a) Time history, (b) Fourier spectra, and (c) phase plane, forf
p
=
0.45.
(a)
(b)
(c)The dependence of the excitationf
p for homoclinic chaos for different values of the damping parameter μ for a frequency range 0
<
ω
<
1.5 has been plotted. The surfaces are shown in Figure 9(a). These surfaces play an important role in the formation of the chaotic region. More precisely, this condition provides a domain on the parameter spaces where the system has transverse homoclinic orbits resulting in possible chaotic behavior. As it can be observed from this figure, as μ increases the threshold f
p obtained by the Melnikov method increases on the frequency 0
<
ω
<
1.5.Homoclinic bifurcation curve for (a) control parameterf
p and (b) control parameter μ.
(a)
(b)Now we consider the effect of the parameterμ for the occurrence of homoclinic bifurcation and chaos in system. Using (17) and choosing μ as control parameter, the condition for transverse intersection of stable and unstable manifolds is given by
(19)
μ
≤
3
2
4
π
ω
c
f
p
b
b
sec
h
(
π
ω
2
b
)
.
In order to have a better understanding, we have plotted in Figure 9(b) the dependence of the parameter μ on the frequency 0
<
ω
<
2 for different values of f
p. In the parameter region below the threshold curve transverse intersection of stable and unstable manifolds occurs, and the onset of chaos is expected. One observation from this figure is that the threshold μ increases when f
p increases. These surfaces provide a domain on the parameter spaces where the system has homoclinic bifurcation resulting in possible chaotic dynamic. More precisely, these results play an important role in the formation of the chaotic regions and could be used for the analysis and design of the gear system parameters.


## 5. Conclusions

In this paper, a nonlinear model of a spur gear system has been proposed to investigate the homoclinic bifurcation and the chaotic behaviors of system. From the Melnikov approach, the threshold values of the control parameter for the occurrence of homoclinic bifurcation and onset of chaos have been predicted. The analytical predictions have been validated with numerical simulations and good agreement is observed. Analyzing and predicting the chaotic behaviors of a gear system are useful. The results provide some idea and guidance for the analysis of the gear system parameters and could be used for designing and developing an optimal gear transmission system.



---

*Source: 809739-2014-05-28.xml*