# Dynamic Analysis of a Fault Planetary Gear System under Nonlinear Parameter Excitation

Authors: Jianchao Han; Yinghui Liu; Lei Liang; Yang Zhao; Huibo Zhang
Journal: Shock and Vibration
 (2021)

---

## Abstract

In order to facilitate lubrication and avoid the gear stuck due to thermal expansion, there needs to be a gap between the tooth profiles. As a strong nonlinear factor, the backlash will affect the motion state of the planetary gear system. When the gear failures occur, the motion state of the system will accordingly change. In this study, the meshing stiffness of the gear pair with tooth tip chipping fault is calculated by combining the analytic geometry method and the potential energy method. Then, a new nonlinear dynamic model including tooth backlash, time-varying mesh stiffness, and manufacturing error is established to study the dynamic response of the system. The equations of motion are derived by the Lagrangian method and solved by the numerical integration method. Taking the excitation frequency and tooth backlash as the variation parameters, respectively, the dynamic characteristics of the system are analyzed by comparing the global bifurcation diagrams between the health system and the fault system, and the path of the system into chaos is revealed. At the same time, the local characteristics of the system are revealed through the phase diagrams and Poincaré maps. The results show that with the variation of excitation frequency and tooth backlash, the fault system presents a more complex motion state. This study can provide the theoretical support for dynamic design and fault diagnosis of planetary gear transmission systems under the environment of gear fault-prone.


---

## Body

## 1. Introduction

Planetary gear systems are widely used in aerospace, agricultural machinery, construction machinery, and other fields because of their high transmission ratio and transmission efficiency [1]. For the purpose that the gears can be fully lubricated and avoid jamming, tooth backlash between the engaged teeth is indispensable. However, the backlash is a strong nonlinear factor for the dynamic characteristics of the gear system. Due to the nonlinear factors, such as backlash, gear fault, and so on, the system exhibits different motion states with the variational excitation frequency. Hence, for better dynamic design and fault diagnosis for planetary gear systems, it is necessary to analyze their dynamic characteristics.In recent decades, many scholars have conducted modeling analyses on the gear transmission system. In 1994, Kahraman [2] proposed a dynamic model of the planetary gear system and studied the load-sharing characteristics. Later, the vibration modes of the planetary gear system [3–5] and the suppression law of the meshing phase to the planetary modal response are analyzed [6]. Then, the load distribution coefficient is deduced based on the model [7–9]. Wang et al. [10] revealed the nonlinear phenomena and evolution mechanisms of bifurcation and chaos via a three-degree-of-freedom torsional vibration model. Then, Shen et al. [11] studied the dynamic characteristics of spur gear systems using the incremental harmonic balance method. In the planetary gear system, gear errors are inevitable, such as installation errors [12], manufacturing errors [13], and geometric errors [14], and the influence on the system is also analyzed [15].Backlash and meshing stiffness are nonlinear factors for the gear system. To analyze the dynamic characteristics of the planetary gear system, the nonlinear dynamic model containing the tooth backlash [16–18] and time-varying meshing stiffness was established [19–22]. Then, Huang et al. [23] and Pan and Vicuña [24] introduced the fractal backlash to the model. However, under high-speed and light-load conditions, gears may appear as tooth backside contact. Therefore, a multistate dynamic model was established by Liu et al. [25], and the variation law of the meshing force under different meshing conditions was analyzed. The influence of bearing clearance [26] and gear surface modification [27] on the dynamic response of the planetary gear system was revealed. Later, a nonlinear dynamic model of the multistage gear transmission was established by Zhao and Ji [28] and Xiang et al. [29], and the nonlinear dynamic characteristics of the system are analyzed at the same time.Due to poor lubrication, impact load, and stress concentration, gear failures often occur and the dynamic response of the system will change accordingly. For early fault detection, Shi et al. [30] analyzed the fault characteristics under variable load. Pan et al. [31] analyzed the frequency components of the vibration signal in the fault and the healthy state. Considering the flexible ring gear and bearing fault, Liu et al. [32] presented a rigid-flexible coupling planetary gear dynamic model. By establishing the meshing stiffness model under spalling fault, Luo et al. [33] and Xiang et al. [34] analyzed the dynamic characteristics of the planetary gear system. Later, Shen et al. [35] proposed a purely torsional model to analyze the dynamic characteristics of planetary gear under wear fault. Yang et al. [36] proposed a nonlinear dynamic model with tooth backlash and bearing clearance and analyzed the vibration response under crack fault [37–39]. Luo et al. [40] established the meshing stiffness model for spalling and pitting failures and compared the dynamic response under different failure types.In the previous literature, the dynamic response of the planetary gear system under crack failure, pitting failure, and spalling was researched. However, there are limited research studies on the nonlinear dynamic characteristics of planetary gear systems with tooth tip chipping fault under nonlinear parameters excitation. Therefore, in order to reveal the dynamic characteristics of the fault system, this study established a nonlinear dynamic model containing tooth backlash, time-varying meshing stiffness, and static transmission error and analyzed the chaos and bifurcation characteristics of the system via choosing the excitation frequency and tooth backlash as variable parameters. The rest of the study is organized as follows. Section2 establishes the planetary gear dynamic model. In Section 3, the dynamic characteristics of the system are analyzed by taking the tooth backlash and rotation speed as control variables. Finally, some conclusions are given in Section 4.


## 2. Dynamic Model of the Planetary Gear System

The planetary gear system has a complex structure that is different from fixed shaft gearboxes. In order to study the nonlinear dynamic behavior, a pure torsional dynamic model of the planetary gear system is proposed, which contains a sun gears, a ring gear r, a carrier c, and N planet gears pi (i = 1, 2, …, N), as shown in Figure 1. In this model, all gears are standard spur gears. The carrier with planet gears fixed on is connected with the input shaft, and the sun gear is attached with the output shaft. The ring gear is fixed. In order to simplify the model, all components are assumed to be rigid.Figure 1
The dynamic model of planetary gear transmission.


### 2.1. System Excitations

#### 2.1.1. Time-Varying Meshing Stiffness of Gears with Tooth Tip Chipping Fault

The potential energy method is a common method for calculating meshing stiffness. In the method, the gear tooth is regarded as a variable cross-section cantilever beam fixed on the tooth root, as shown in Figure2. The total meshing stiffness can be composed of bending stiffness kb, axial compression stiffness ka, shear stiffness, Hertz contact stiffness kh, and fillet foundation stiffness kf. When in the gear occurs tooth tip chipping failure, a part of the material will fall off the tooth tip, which will affect the cross-sectional area and the moment of inertia of the tooth faulty part, as shown in Figure 3. Thereby, to estimate the meshing stiffness of gear pairs with tooth tip chipping, an accurate stiffness model is established.Figure 2
Cantilever beam of spur gear with tooth tip chipping.Figure 3
Diagram of the tooth tip chipping part for the tooth.
(a)(b)In order to simplify the model, the fracture surface is simplified as a plane, as shown in Figure3. The intersection line between the plane and the tooth profile surface is the curve L3, and the intersection line with the gear end surface denotes the straight L2. In the three-dimensional coordinate system, the involute tooth profile equation Ω1 can be expressed as(1)xyz=cosα2−sinα20sinα2cosα20000uvw,where u and v represent the coordinate of the involute in the UOV coordinate system, which can be described as(2)u=Rbsinθx+αx−θx+αxcosθx+αx,v=Rbcosθx+αx+θx+αxsinθx+αx.In the three-dimensional coordinate system, the fracture surface equationΩ2 can be written as(3)z=x−xaya−yb−xa−xby−yaxa−xbya−yc+ya−ybxc−xazc.The equation of the curveL2 can be deduced by simultaneous equations of (1) and (3). The equation of line L3 in Figure 3(a) can be expressed as(4)x=Rb/cosαasinβa−xbRb/cosαacosβa−m/2ze+2ha∗cosβty−Rbcosαacosβa+Rbcosαasinβa,βa=π2ze−tanαa−αa−tanα0−α0,βt=π2ze−tanαt−αt−tanα0−α0.As tooth tip chipping failure occurs, a part of the material will fall off the tooth tip, so the cross-sectional area and moment of inertia of the faulty part will vary accordingly, which can be deduced by the following equations:(5)Ax′=2hxL−12d1d2,Ix′=Ix−d2d1336+12d1d2hx−13d12−1/2d1d2hx−1/3d122hxL−1/2d1d2,d1=yxL1−yxL3,d2=zxL2.Thereby, the meshing stiffness can be obtained by the potential energy method [41]:(6)1ka=∫0dsinα12EAxdx,1kb=∫0dd−xcosα1−hsinα12EIxdx,1ks=∫0d1.2cosα122GAxdx,kh=πEL−d241−ν2,1kf=cos2α1ELL∗ufSf2+M∗ufSf+P∗1+Q∗tan2α1,where the parameters of L∗, M∗, P∗, Q∗, uf, and Sf are given in the study by Ma et al. [42].The total meshing stiffness of external-external mesh can be expressed as(7)kspn=∑i=1211/kh,i+1/kb1,i+1/ks1,i+1/ka1,i+1/kf1,i+1/kb2,i+1/ks2,i+1/ka2,i+1/kf2,i.In this model, the internal-external mesh is assumed to be in a healthy state, and the total meshing stiffness can be deduced in the same way:(8)krpn=∑i=1n11/kh,i+1/kb1,i+1/ks1,i+1/ka1,i+1/kb2,i+1/ks2,i+1/ka2,i+1/kf2,i,where n represents the number of tooth pairs engaged at the same time. The subscript 1 indicates the driving gear, and the subscript 2 represents the driven gear.


#### 2.1.2. Manufacturing Error and Damping

In this model, the carrier is connected to the input shaft, and the meshing frequency can be expressed as(9)ωm=ωczr.Assuming that the gear teeth are exactly the same, the static transmission error between the engaged gear teeth can be described as Fourier series, where the meshing frequency is the fundamental frequency. To simplify the model, merely, the fundamental frequency is taken into account, which can be expressed as [29](10a)espn=easpnsinωmt+ϕspn,(10b)erpn=earpnsinωmt+ϕrpn.The mesh damping can be expressed as(11a)cspn=2ξspnkmspn1/ms+1/mpn,(11b)crpn=2ξrpnkmrpn1/mr+1/mpn,where kmspn and kmrpn are the arithmetic mean of kspn and krpn, Ij (j=s,ce,r,pn) is the moment of inertia, and mj (j=s,c,r,pn) is the equivalent mass that can be calculated by the following equations:(12)Ice=Ic+Nmpnrc2,ms=Isrs2,mc=Icerc2,mpn=Ipnrpn2.


#### 2.1.3. Gear Backlash

In order to facilitate gear lubrication and prevent the gear from jamming during the operation, there is a gap between the engaged gear teeth. The backlash can be described by a piecewise function:(13a)fδspn=δspn−bspn,δspn>bspn,0,δspn≤bspn,δspn+bspn,δspn<−bspn,(13b)fδrpn=δrpn−brpn,δrpn>brpn,0,δrpn≤brpn,δrpn+brpn,δrpn<−brpn.


### 2.2. Dynamic Differential Equations

According to the Lagrangian equation, the differential equations of torsional vibration under the excitations can be obtained:(14)Isθ¨s+∑n=1NrsFspn=Tout,Ipnθ¨pn−rpnFspn+rpnFrpn=0,Iceθ¨c−∑n=1NrcFspn−∑n=1NrcFrpn=−Tin,where Tout and Tin represent the output torque and the input torque, respectively. Fspn and Frpn are adopted to describe the meshing force which can be expressed as [20](15a)Fspn=cspnδ˙spn+kspnfδspn,(15b)Frpn=crpnδ˙rpn+krpnfδrpn.In order to reduce the dimension of the differential equations, generalized coordinates are introduced as follows:(16a)δspn=rsθs−rpnθpn−rcθc−espn,(16b)δrpn=rpnθpn−rrθr−rcθc−erpn.Due to the large numerical difference between the various nonlinear parameters, the solving progress is difficult to conduct. In order to facilitate the solution and improve the calculation efficiency, the dimensionless processing is carried out. By introducing the time scaleωn and displacement scale bc, the dimensionless parameters can be obtained. The dimensionless time displacement, velocity, and acceleration can be expressed as(17)τ=ωn∗t,δ¯spn=δspnbc,δ¯˙spn=δspnbcωn,δ¯¨spn=δ¯¨spnbcωn2,δ¯rpn=δrpnbc,δ¯˙rpn=δrpnbcωn,δ¯¨rpn=δ¯¨rpnbcωn2,ωn=kmrpn1/ms+1/mc.The dimensionless mesh stiffness, static transmission error, and backlash are expressed as(18a)e¯spn=easpnbcsinωmωnt+ψspn,(18b)e¯rpn=earpnbcsinωmωnt+ψspn,(19a)fδ¯spn=δ¯spn−bspnbc,δ¯spn>bspnbc,0,δ¯spn≤bspnbc,δ¯spn+bspnbc,δ¯spn<−bspnbc,(19b)fδ¯rpn=δ¯rpn−brpnbc,δ¯rpn>brpnbc,0,δ¯rpn≤brpnbc,δ¯rpn+brpnbc,δ¯rpn<−brpnbc.By substituting equations (15a)–(19b) into equation (14), the differential equation of dimensionless parameters can be obtained as follows:(20)δ¯¨spn+1ms+1mc1ωn∑n=1Ncspnδ¯˙spn+1ms+1mc1ωn2∑n=1Nk¯spnfδ¯spn+1mpn1ωncspnδ¯˙spn−1mpn1ωncrpnδ¯˙rpn−1mc1ωn∑n=1Ncrpnδ¯˙rpn+1mc1ωn2∑n=1Nk¯rpnfδ¯spn=ToutrsIsbcωn2+TinrcIcbcωn2+easpnbcωmωn2sinωmωnt+ψspn,δ¯¨rpn−1mpn1ωncspnδ¯˙spn−1mpn1ωn2k¯spnfδ¯spn+1mpn1ωncrpnδ¯˙rpn+1mc1ωn2k¯rpnfδ¯˙rpn+1mc1ωn∑n=1Ncspnδ¯˙spn+1mc1ωn2∑n=1Nk¯spnfδ¯spn+1mc1ωncrpnδ¯˙rpn+1mc1ωn2∑n=1Nk¯spnfδ¯spn=TinrcIcbcωn2−earpnbcωmωn2sinωmωnt+ψspn.


## 2.1. System Excitations

### 2.1.1. Time-Varying Meshing Stiffness of Gears with Tooth Tip Chipping Fault

The potential energy method is a common method for calculating meshing stiffness. In the method, the gear tooth is regarded as a variable cross-section cantilever beam fixed on the tooth root, as shown in Figure2. The total meshing stiffness can be composed of bending stiffness kb, axial compression stiffness ka, shear stiffness, Hertz contact stiffness kh, and fillet foundation stiffness kf. When in the gear occurs tooth tip chipping failure, a part of the material will fall off the tooth tip, which will affect the cross-sectional area and the moment of inertia of the tooth faulty part, as shown in Figure 3. Thereby, to estimate the meshing stiffness of gear pairs with tooth tip chipping, an accurate stiffness model is established.Figure 2
Cantilever beam of spur gear with tooth tip chipping.Figure 3
Diagram of the tooth tip chipping part for the tooth.
(a)(b)In order to simplify the model, the fracture surface is simplified as a plane, as shown in Figure3. The intersection line between the plane and the tooth profile surface is the curve L3, and the intersection line with the gear end surface denotes the straight L2. In the three-dimensional coordinate system, the involute tooth profile equation Ω1 can be expressed as(1)xyz=cosα2−sinα20sinα2cosα20000uvw,where u and v represent the coordinate of the involute in the UOV coordinate system, which can be described as(2)u=Rbsinθx+αx−θx+αxcosθx+αx,v=Rbcosθx+αx+θx+αxsinθx+αx.In the three-dimensional coordinate system, the fracture surface equationΩ2 can be written as(3)z=x−xaya−yb−xa−xby−yaxa−xbya−yc+ya−ybxc−xazc.The equation of the curveL2 can be deduced by simultaneous equations of (1) and (3). The equation of line L3 in Figure 3(a) can be expressed as(4)x=Rb/cosαasinβa−xbRb/cosαacosβa−m/2ze+2ha∗cosβty−Rbcosαacosβa+Rbcosαasinβa,βa=π2ze−tanαa−αa−tanα0−α0,βt=π2ze−tanαt−αt−tanα0−α0.As tooth tip chipping failure occurs, a part of the material will fall off the tooth tip, so the cross-sectional area and moment of inertia of the faulty part will vary accordingly, which can be deduced by the following equations:(5)Ax′=2hxL−12d1d2,Ix′=Ix−d2d1336+12d1d2hx−13d12−1/2d1d2hx−1/3d122hxL−1/2d1d2,d1=yxL1−yxL3,d2=zxL2.Thereby, the meshing stiffness can be obtained by the potential energy method [41]:(6)1ka=∫0dsinα12EAxdx,1kb=∫0dd−xcosα1−hsinα12EIxdx,1ks=∫0d1.2cosα122GAxdx,kh=πEL−d241−ν2,1kf=cos2α1ELL∗ufSf2+M∗ufSf+P∗1+Q∗tan2α1,where the parameters of L∗, M∗, P∗, Q∗, uf, and Sf are given in the study by Ma et al. [42].The total meshing stiffness of external-external mesh can be expressed as(7)kspn=∑i=1211/kh,i+1/kb1,i+1/ks1,i+1/ka1,i+1/kf1,i+1/kb2,i+1/ks2,i+1/ka2,i+1/kf2,i.In this model, the internal-external mesh is assumed to be in a healthy state, and the total meshing stiffness can be deduced in the same way:(8)krpn=∑i=1n11/kh,i+1/kb1,i+1/ks1,i+1/ka1,i+1/kb2,i+1/ks2,i+1/ka2,i+1/kf2,i,where n represents the number of tooth pairs engaged at the same time. The subscript 1 indicates the driving gear, and the subscript 2 represents the driven gear.


### 2.1.2. Manufacturing Error and Damping

In this model, the carrier is connected to the input shaft, and the meshing frequency can be expressed as(9)ωm=ωczr.Assuming that the gear teeth are exactly the same, the static transmission error between the engaged gear teeth can be described as Fourier series, where the meshing frequency is the fundamental frequency. To simplify the model, merely, the fundamental frequency is taken into account, which can be expressed as [29](10a)espn=easpnsinωmt+ϕspn,(10b)erpn=earpnsinωmt+ϕrpn.The mesh damping can be expressed as(11a)cspn=2ξspnkmspn1/ms+1/mpn,(11b)crpn=2ξrpnkmrpn1/mr+1/mpn,where kmspn and kmrpn are the arithmetic mean of kspn and krpn, Ij (j=s,ce,r,pn) is the moment of inertia, and mj (j=s,c,r,pn) is the equivalent mass that can be calculated by the following equations:(12)Ice=Ic+Nmpnrc2,ms=Isrs2,mc=Icerc2,mpn=Ipnrpn2.


### 2.1.3. Gear Backlash

In order to facilitate gear lubrication and prevent the gear from jamming during the operation, there is a gap between the engaged gear teeth. The backlash can be described by a piecewise function:(13a)fδspn=δspn−bspn,δspn>bspn,0,δspn≤bspn,δspn+bspn,δspn<−bspn,(13b)fδrpn=δrpn−brpn,δrpn>brpn,0,δrpn≤brpn,δrpn+brpn,δrpn<−brpn.


## 2.1.1. Time-Varying Meshing Stiffness of Gears with Tooth Tip Chipping Fault

The potential energy method is a common method for calculating meshing stiffness. In the method, the gear tooth is regarded as a variable cross-section cantilever beam fixed on the tooth root, as shown in Figure2. The total meshing stiffness can be composed of bending stiffness kb, axial compression stiffness ka, shear stiffness, Hertz contact stiffness kh, and fillet foundation stiffness kf. When in the gear occurs tooth tip chipping failure, a part of the material will fall off the tooth tip, which will affect the cross-sectional area and the moment of inertia of the tooth faulty part, as shown in Figure 3. Thereby, to estimate the meshing stiffness of gear pairs with tooth tip chipping, an accurate stiffness model is established.Figure 2
Cantilever beam of spur gear with tooth tip chipping.Figure 3
Diagram of the tooth tip chipping part for the tooth.
(a)(b)In order to simplify the model, the fracture surface is simplified as a plane, as shown in Figure3. The intersection line between the plane and the tooth profile surface is the curve L3, and the intersection line with the gear end surface denotes the straight L2. In the three-dimensional coordinate system, the involute tooth profile equation Ω1 can be expressed as(1)xyz=cosα2−sinα20sinα2cosα20000uvw,where u and v represent the coordinate of the involute in the UOV coordinate system, which can be described as(2)u=Rbsinθx+αx−θx+αxcosθx+αx,v=Rbcosθx+αx+θx+αxsinθx+αx.In the three-dimensional coordinate system, the fracture surface equationΩ2 can be written as(3)z=x−xaya−yb−xa−xby−yaxa−xbya−yc+ya−ybxc−xazc.The equation of the curveL2 can be deduced by simultaneous equations of (1) and (3). The equation of line L3 in Figure 3(a) can be expressed as(4)x=Rb/cosαasinβa−xbRb/cosαacosβa−m/2ze+2ha∗cosβty−Rbcosαacosβa+Rbcosαasinβa,βa=π2ze−tanαa−αa−tanα0−α0,βt=π2ze−tanαt−αt−tanα0−α0.As tooth tip chipping failure occurs, a part of the material will fall off the tooth tip, so the cross-sectional area and moment of inertia of the faulty part will vary accordingly, which can be deduced by the following equations:(5)Ax′=2hxL−12d1d2,Ix′=Ix−d2d1336+12d1d2hx−13d12−1/2d1d2hx−1/3d122hxL−1/2d1d2,d1=yxL1−yxL3,d2=zxL2.Thereby, the meshing stiffness can be obtained by the potential energy method [41]:(6)1ka=∫0dsinα12EAxdx,1kb=∫0dd−xcosα1−hsinα12EIxdx,1ks=∫0d1.2cosα122GAxdx,kh=πEL−d241−ν2,1kf=cos2α1ELL∗ufSf2+M∗ufSf+P∗1+Q∗tan2α1,where the parameters of L∗, M∗, P∗, Q∗, uf, and Sf are given in the study by Ma et al. [42].The total meshing stiffness of external-external mesh can be expressed as(7)kspn=∑i=1211/kh,i+1/kb1,i+1/ks1,i+1/ka1,i+1/kf1,i+1/kb2,i+1/ks2,i+1/ka2,i+1/kf2,i.In this model, the internal-external mesh is assumed to be in a healthy state, and the total meshing stiffness can be deduced in the same way:(8)krpn=∑i=1n11/kh,i+1/kb1,i+1/ks1,i+1/ka1,i+1/kb2,i+1/ks2,i+1/ka2,i+1/kf2,i,where n represents the number of tooth pairs engaged at the same time. The subscript 1 indicates the driving gear, and the subscript 2 represents the driven gear.


## 2.1.2. Manufacturing Error and Damping

In this model, the carrier is connected to the input shaft, and the meshing frequency can be expressed as(9)ωm=ωczr.Assuming that the gear teeth are exactly the same, the static transmission error between the engaged gear teeth can be described as Fourier series, where the meshing frequency is the fundamental frequency. To simplify the model, merely, the fundamental frequency is taken into account, which can be expressed as [29](10a)espn=easpnsinωmt+ϕspn,(10b)erpn=earpnsinωmt+ϕrpn.The mesh damping can be expressed as(11a)cspn=2ξspnkmspn1/ms+1/mpn,(11b)crpn=2ξrpnkmrpn1/mr+1/mpn,where kmspn and kmrpn are the arithmetic mean of kspn and krpn, Ij (j=s,ce,r,pn) is the moment of inertia, and mj (j=s,c,r,pn) is the equivalent mass that can be calculated by the following equations:(12)Ice=Ic+Nmpnrc2,ms=Isrs2,mc=Icerc2,mpn=Ipnrpn2.


## 2.1.3. Gear Backlash

In order to facilitate gear lubrication and prevent the gear from jamming during the operation, there is a gap between the engaged gear teeth. The backlash can be described by a piecewise function:(13a)fδspn=δspn−bspn,δspn>bspn,0,δspn≤bspn,δspn+bspn,δspn<−bspn,(13b)fδrpn=δrpn−brpn,δrpn>brpn,0,δrpn≤brpn,δrpn+brpn,δrpn<−brpn.


## 2.2. Dynamic Differential Equations

According to the Lagrangian equation, the differential equations of torsional vibration under the excitations can be obtained:(14)Isθ¨s+∑n=1NrsFspn=Tout,Ipnθ¨pn−rpnFspn+rpnFrpn=0,Iceθ¨c−∑n=1NrcFspn−∑n=1NrcFrpn=−Tin,where Tout and Tin represent the output torque and the input torque, respectively. Fspn and Frpn are adopted to describe the meshing force which can be expressed as [20](15a)Fspn=cspnδ˙spn+kspnfδspn,(15b)Frpn=crpnδ˙rpn+krpnfδrpn.In order to reduce the dimension of the differential equations, generalized coordinates are introduced as follows:(16a)δspn=rsθs−rpnθpn−rcθc−espn,(16b)δrpn=rpnθpn−rrθr−rcθc−erpn.Due to the large numerical difference between the various nonlinear parameters, the solving progress is difficult to conduct. In order to facilitate the solution and improve the calculation efficiency, the dimensionless processing is carried out. By introducing the time scaleωn and displacement scale bc, the dimensionless parameters can be obtained. The dimensionless time displacement, velocity, and acceleration can be expressed as(17)τ=ωn∗t,δ¯spn=δspnbc,δ¯˙spn=δspnbcωn,δ¯¨spn=δ¯¨spnbcωn2,δ¯rpn=δrpnbc,δ¯˙rpn=δrpnbcωn,δ¯¨rpn=δ¯¨rpnbcωn2,ωn=kmrpn1/ms+1/mc.The dimensionless mesh stiffness, static transmission error, and backlash are expressed as(18a)e¯spn=easpnbcsinωmωnt+ψspn,(18b)e¯rpn=earpnbcsinωmωnt+ψspn,(19a)fδ¯spn=δ¯spn−bspnbc,δ¯spn>bspnbc,0,δ¯spn≤bspnbc,δ¯spn+bspnbc,δ¯spn<−bspnbc,(19b)fδ¯rpn=δ¯rpn−brpnbc,δ¯rpn>brpnbc,0,δ¯rpn≤brpnbc,δ¯rpn+brpnbc,δ¯rpn<−brpnbc.By substituting equations (15a)–(19b) into equation (14), the differential equation of dimensionless parameters can be obtained as follows:(20)δ¯¨spn+1ms+1mc1ωn∑n=1Ncspnδ¯˙spn+1ms+1mc1ωn2∑n=1Nk¯spnfδ¯spn+1mpn1ωncspnδ¯˙spn−1mpn1ωncrpnδ¯˙rpn−1mc1ωn∑n=1Ncrpnδ¯˙rpn+1mc1ωn2∑n=1Nk¯rpnfδ¯spn=ToutrsIsbcωn2+TinrcIcbcωn2+easpnbcωmωn2sinωmωnt+ψspn,δ¯¨rpn−1mpn1ωncspnδ¯˙spn−1mpn1ωn2k¯spnfδ¯spn+1mpn1ωncrpnδ¯˙rpn+1mc1ωn2k¯rpnfδ¯˙rpn+1mc1ωn∑n=1Ncspnδ¯˙spn+1mc1ωn2∑n=1Nk¯spnfδ¯spn+1mc1ωncrpnδ¯˙rpn+1mc1ωn2∑n=1Nk¯spnfδ¯spn=TinrcIcbcωn2−earpnbcωmωn2sinωmωnt+ψspn.


## 3. Numerical Simulation and Results Analysis

Since the relative displacement between the engaged teeth has the same law, the dimensionless relative displacement of the engaged teeth between the sun gear and the planet gear is taken as the object to analyze, and the fourth-order Runge–Kutta method is used to solve the dynamics differential equations. The basic parameters of each component in the model are given in Table1, and the meshing parameters are given in Table 2. The displacement scale bc is 1e − 5 m. The input torque is 100 Nm, and the transmission ratio is 4. The αa is assigned a value of 25°, yb is 0 mm, and zc is 20 mm. The meshing stiffness of external-external mesh and internal-external mesh is shown in Figures 4 and 5. To avoid the influence of transient response, the first 500 response cycles of the system are omitted.Table 1
Parameters of the dynamic model.
Parameter nameSun gearRing gearPlanet gearThe carrierNumber of parts1131Number of teeth195719—Gear module (mm)3.23.23.2—Pressure angle (°)202020—Inertia moment (kg·m^2)3.85E − 042.4E − 023.85E − 044.93E − 03Table 2
Meshing parameters.
Parameter nameSun-planetRing-planetDamping ratio0.070.07Error amplitude (m)1e − 51e − 5Figure 4
Meshing stiffness of external-external mesh.Figure 5
Meshing stiffness of internal-external mesh.


### 3.1. Bifurcation and Chaos of Fault and Health Systems with Excitation Frequency

The motion state of the system will convert with the variation of the excitation frequency. Therefore, the dimensionless excitation frequency Ω is chosen as the variable parameter, and the dimensionless backlashb is assigned the value of 2; the bifurcation diagram of the system in both health and fault conditions is shown in Figures 6 and 7. The largest Lyapunov exponent diagram of the health system is shown in Figure 8. As can be seen from figures, both the health and fault systems have rich chaotic and bifurcation characteristics. For the health system, when the excitation frequency Ω is between 0 and 0.69, the system presents a single-period motion and the corresponding largest Lyapunov exponent is less than zero. When it is between 0.69 and 1.28, the system undergoes a single-period motion to evolve into chaotic motion with small amplitude. When the excitation frequency changes in the range of 1.28–1.48, the system basically keeps cyclical movement, in which the largest Lyapunov exponent is less than zero. Then, the system enters the chaotic motion eventually. For the fault system, it is similar to the state of motion experienced by the health system, but in the excitation frequency range corresponding to the periodic motion of the health system, the fault system shows obvious multiple periodic motions, as shown in Figure 9. Under different tooth backlashes, the bifurcation diagram of the system is displayed in Figure 10. It can be seen that the amplitude of chaotic motion augments with the increase of the initial backlash. The system enters the chaotic motion earlier as the speed increases if the backlash is designed large.Figure 6
Bifurcation diagram of the health system with excitation frequency variation.Figure 7
Bifurcation diagram of the fault system with excitation frequency variation.Figure 8
Largest Lyapunov exponent diagram of the health system with frequency variation.Figure 9
Bifurcation comparison diagram with frequency of the fault system and health system.Figure 10
Bifurcation comparison diagram with frequency of the fault system and health system.In order to further analyze the local characteristics of the health system and fault system, phase diagrams and Poincaré maps are used to describe the local characteristics of the system. Figures11–15 are the phase diagram and Poincaré maps of the two systems at different excitation frequencies. It can be seen from Figure 11 that the health system exhibits a single-period motion at Ω = 0.5, while the fault system exhibits a multiple period motion. When Ω equals to 1, as shown in Figure 12, both of the two systems are in chaotic motion. When Ω is 1.34, as shown in Figure 13, the health system exhibits single-period motion, while the fault system presents multiple periodic motions. When Ω is 4, as shown in Figure 14, the health system exhibits quasiperiodic motion, while the faulty system appears as small amplitude chaotic motion. As Ω continues to increase, the health system exhibits the single-periodic motion at Ω = 1.47, while the faulty system occurs in chaotic motion, as shown in Figure 15. Hence, when the health system exhibits periodic motion, the fault system exhibits the multiperiod or pseudoperiodic motion at the corresponding excitation frequency.Figure 11
(a) Health system atΩ=0.5. (b) Fault system at Ω=0.5.
(a)(b)Figure 12
(a) Health system atΩ=1. (b) Fault system at Ω=1.
(a)(b)Figure 13
(a) Health system atΩ=1.34. (b) Fault system at Ω=1.34.
(a)(b)Figure 14
(a) Health system atΩ=1.4. (b) Fault system at Ω=1.4.
(a)(b)Figure 15
(a) Health system atΩ=1.47. (b) Fault system at Ω=1.47.
(a)(b)


### 3.2. Bifurcation and Chaos of Fault and Health Systems with Different Initial Backlashes

In order to facilitate lubrication, there needs to be a gap between the engaged gear teeth and the value of the initial backlash that will affect the dynamic characteristics of the gear system. Therefore, the dimensionless backlashb is chosen as the variable parameter to study the dynamic characteristics of the system with the dimensionless excitation frequency Ω to be 1.6. When the backlash varies between 0 and 3, the bifurcation diagrams of the health system and the fault system are shown in Figures 16 and 17, respectively. The largest Lyapunov exponent diagram of the health system is shown in Figure 18. It can be seen from the figures that the amplitude of the system motion augments with the increase of the backlash. For the health system, when the backlash b changes between 0 and 1.095, the system presents a single-period motion and the largest Lyapunov exponent is less than zero. However, the system motion state occurs jump variation at b = 1.095. Then, the system enters a small amplitude chaotic motion at the interval of 1.095–1.29 corresponding to which the largest Lyapunov exponent is larger than zero. Afterwards, the system reenters the single-period motion when b varies in 1.29–1.53. As b continues to increase, and the system eventually enters chaotic motion. For the fault system, it is similar to the state of motion experienced by the health system, but in the backlash interval corresponding to the periodic motion state of the health system, the fault system shows obvious small amplitude chaotic motion. Figure 19 shows the bifurcation characteristics of the two health systems with backlash. It can be seen from the figure that the fault system has more complex nonlinear characteristics and enters the chaotic motion earlier. After changing the speed, the bifurcation diagram of the fault system with the backlash is shown in Figure 20. When the speed increases, the amplitude of the chaotic motion will increase accordingly.Figure 16
Bifurcation diagram of the health system with backlash variation.Figure 17
Bifurcation diagram of the fault system with backlash variation.Figure 18
Largest Lyapunov exponent diagram of the health system with backlash variation.Figure 19
Bifurcation comparison diagram with backlash of the fault system and health system.Figure 20
Bifurcation comparison diagram with backlash of the fault system and health system.In order to further analyze the local characteristics of the health system and fault system and reveal the path of the system into chaos, phase diagrams and Poincaré maps are used to describe the local characteristics of the system. From Figures21–24, it can be seen that the health system exhibits a single-period motion at b = 0.5, 1, and 1.4 and a small amplitude chaotic motion at 1.2. While, the fault system directly enters the chaotic motion.Figure 21
(a) Health system atb=0.5. (b) Fault system at b=0.5.
(a)(b)Figure 22
(a) Health system atb=1. (b) Fault system at b=1.
(a)(b)Figure 23
(a) Health system atb=1.2. (b) Fault system at b=1.2.
(a)(b)Figure 24
(a) Health system atb=1.4. (b) Fault system at b=1.4.
(a)(b)


## 3.1. Bifurcation and Chaos of Fault and Health Systems with Excitation Frequency

The motion state of the system will convert with the variation of the excitation frequency. Therefore, the dimensionless excitation frequency Ω is chosen as the variable parameter, and the dimensionless backlashb is assigned the value of 2; the bifurcation diagram of the system in both health and fault conditions is shown in Figures 6 and 7. The largest Lyapunov exponent diagram of the health system is shown in Figure 8. As can be seen from figures, both the health and fault systems have rich chaotic and bifurcation characteristics. For the health system, when the excitation frequency Ω is between 0 and 0.69, the system presents a single-period motion and the corresponding largest Lyapunov exponent is less than zero. When it is between 0.69 and 1.28, the system undergoes a single-period motion to evolve into chaotic motion with small amplitude. When the excitation frequency changes in the range of 1.28–1.48, the system basically keeps cyclical movement, in which the largest Lyapunov exponent is less than zero. Then, the system enters the chaotic motion eventually. For the fault system, it is similar to the state of motion experienced by the health system, but in the excitation frequency range corresponding to the periodic motion of the health system, the fault system shows obvious multiple periodic motions, as shown in Figure 9. Under different tooth backlashes, the bifurcation diagram of the system is displayed in Figure 10. It can be seen that the amplitude of chaotic motion augments with the increase of the initial backlash. The system enters the chaotic motion earlier as the speed increases if the backlash is designed large.Figure 6
Bifurcation diagram of the health system with excitation frequency variation.Figure 7
Bifurcation diagram of the fault system with excitation frequency variation.Figure 8
Largest Lyapunov exponent diagram of the health system with frequency variation.Figure 9
Bifurcation comparison diagram with frequency of the fault system and health system.Figure 10
Bifurcation comparison diagram with frequency of the fault system and health system.In order to further analyze the local characteristics of the health system and fault system, phase diagrams and Poincaré maps are used to describe the local characteristics of the system. Figures11–15 are the phase diagram and Poincaré maps of the two systems at different excitation frequencies. It can be seen from Figure 11 that the health system exhibits a single-period motion at Ω = 0.5, while the fault system exhibits a multiple period motion. When Ω equals to 1, as shown in Figure 12, both of the two systems are in chaotic motion. When Ω is 1.34, as shown in Figure 13, the health system exhibits single-period motion, while the fault system presents multiple periodic motions. When Ω is 4, as shown in Figure 14, the health system exhibits quasiperiodic motion, while the faulty system appears as small amplitude chaotic motion. As Ω continues to increase, the health system exhibits the single-periodic motion at Ω = 1.47, while the faulty system occurs in chaotic motion, as shown in Figure 15. Hence, when the health system exhibits periodic motion, the fault system exhibits the multiperiod or pseudoperiodic motion at the corresponding excitation frequency.Figure 11
(a) Health system atΩ=0.5. (b) Fault system at Ω=0.5.
(a)(b)Figure 12
(a) Health system atΩ=1. (b) Fault system at Ω=1.
(a)(b)Figure 13
(a) Health system atΩ=1.34. (b) Fault system at Ω=1.34.
(a)(b)Figure 14
(a) Health system atΩ=1.4. (b) Fault system at Ω=1.4.
(a)(b)Figure 15
(a) Health system atΩ=1.47. (b) Fault system at Ω=1.47.
(a)(b)


## 3.2. Bifurcation and Chaos of Fault and Health Systems with Different Initial Backlashes

In order to facilitate lubrication, there needs to be a gap between the engaged gear teeth and the value of the initial backlash that will affect the dynamic characteristics of the gear system. Therefore, the dimensionless backlashb is chosen as the variable parameter to study the dynamic characteristics of the system with the dimensionless excitation frequency Ω to be 1.6. When the backlash varies between 0 and 3, the bifurcation diagrams of the health system and the fault system are shown in Figures 16 and 17, respectively. The largest Lyapunov exponent diagram of the health system is shown in Figure 18. It can be seen from the figures that the amplitude of the system motion augments with the increase of the backlash. For the health system, when the backlash b changes between 0 and 1.095, the system presents a single-period motion and the largest Lyapunov exponent is less than zero. However, the system motion state occurs jump variation at b = 1.095. Then, the system enters a small amplitude chaotic motion at the interval of 1.095–1.29 corresponding to which the largest Lyapunov exponent is larger than zero. Afterwards, the system reenters the single-period motion when b varies in 1.29–1.53. As b continues to increase, and the system eventually enters chaotic motion. For the fault system, it is similar to the state of motion experienced by the health system, but in the backlash interval corresponding to the periodic motion state of the health system, the fault system shows obvious small amplitude chaotic motion. Figure 19 shows the bifurcation characteristics of the two health systems with backlash. It can be seen from the figure that the fault system has more complex nonlinear characteristics and enters the chaotic motion earlier. After changing the speed, the bifurcation diagram of the fault system with the backlash is shown in Figure 20. When the speed increases, the amplitude of the chaotic motion will increase accordingly.Figure 16
Bifurcation diagram of the health system with backlash variation.Figure 17
Bifurcation diagram of the fault system with backlash variation.Figure 18
Largest Lyapunov exponent diagram of the health system with backlash variation.Figure 19
Bifurcation comparison diagram with backlash of the fault system and health system.Figure 20
Bifurcation comparison diagram with backlash of the fault system and health system.In order to further analyze the local characteristics of the health system and fault system and reveal the path of the system into chaos, phase diagrams and Poincaré maps are used to describe the local characteristics of the system. From Figures21–24, it can be seen that the health system exhibits a single-period motion at b = 0.5, 1, and 1.4 and a small amplitude chaotic motion at 1.2. While, the fault system directly enters the chaotic motion.Figure 21
(a) Health system atb=0.5. (b) Fault system at b=0.5.
(a)(b)Figure 22
(a) Health system atb=1. (b) Fault system at b=1.
(a)(b)Figure 23
(a) Health system atb=1.2. (b) Fault system at b=1.2.
(a)(b)Figure 24
(a) Health system atb=1.4. (b) Fault system at b=1.4.
(a)(b)


## 4. Conclusions

This study proposed a meshing stiffness model of gear pair with tooth tip chipping fault by combining the analytic geometry method and potential energy method. Then, a new nonlinear dynamic model of the planetary gear system under health and fault conditions is established considering the time-varying mesh stiffness, tooth backlash, and static transmission error. The model is derived by the Lagrangian method and solved by the numerical integration method. Taking the dynamic transmission error between the sun gear and planet gear as the research object, through the global bifurcation diagram, the variations of the two system motion states with the excitation frequency and tooth backlash are analyzed, and the local characteristics of the systems are analyzed via the phase diagrams and the Poincaré maps. Some conclusions can be obtained as follows:(1)
With the variation of excitation frequency, both the health system and fault system show complex bifurcation and chaotic characteristics, and the chaotic motion is mixed with periodic windows;(2)
The fault system has more complex nonlinear characteristics and enters the chaotic motion earlier. The fault system shows multiple periodic motions and a small amplitude chaotic motion in the interval of excitation frequency and backlash, in which the health system is in the periodic motion state;(3)
The amplitude of the chaotic motion augments with the increase of the initial backlash. The system enters the chaotic motion earlier as the speed increases if the backlash is designed large;(4)
When the initial backlash of the design is determined, the amplitude of the chaotic motion will increase with the speed



---

*Source: 1787525-2021-07-14.xml*