Dynamic Modeling and Nonlinear Analysis of a Spur Gear System Considering a Nonuniformly Distributed Meshing Force

Abstract

In previous studies, the meshing force of a gear system is usually treated as being uniformly distributed for the convenience of analysis. In practical applications, however, it is nonuniformly distributed along the line of action due to machining errors, assembly errors, misalignment errors, etc. When a nonuniformly distributed meshing force is coupled with the shaft deformation, dynamic center distance, and time-varying meshing stiffness, the transmission performance of the gear system will be seriously degraded. Therefore, a nonuniformly distributed meshing force cannot be ignored when considering the gear systems used in complicated working conditions. In this study, the gear’s nonuniformly distributed meshing force is analyzed. Then, an 18 degrees-of-freedom bending-torsion-swing-coupled dynamic model of a pair of involute spur gears is put forward. Through this model, the coupling relationship between the nonuniformly distributed meshing force, shaft bending deformation, and dynamic center distance is accurately described. The influence of meshing frequency, stiffness excitation, damping, and error excitation on the nonlinear dynamic characteristics of the gear system was researched through bifurcation diagrams, phase diagrams, Poincaré maps, and time-domain diagrams. Various complicated nonlinear dynamic behaviors, such as quasiperiodic motion, bifurcation, chaotic motion, and chaotic banding, are revealed. Reasonable parameter ranges that guarantee the gear system is in a stable motion were extracted. By evading complicated nonlinear dynamic behavior, the transmission performance of a gear system was improved. This research will contribute to reducing the vibration and noise of gear systems.

1. Introduction

Gear systems have been widely used in many industrial applications due to their advantages of having accurate transmission ratios, compact dimensions, and high efficiency. With increasing requirements for gear performance in cutting-edge fields, many nonlinear dynamic characteristics that affect transmission performance must be considered in gear design. Therefore, it is important to establish a relatively accurate dynamic model and investigate the complicated nonlinear behavior of the gear system subjected to various excitations.

In recent years, numerous well-established dynamic models have been put forward. Ren [1] analyzed the dynamic behavior of gears caused by translation-torsion coupled vibration. Han et al. [2] deduced the dynamic equations of the spur gear from an energy point of view. In addition to the torsional vibration, the pitch and yaw vibration were also considered in their model. Yang [3] studied the vibration response of a multistage gear system under deterministic and random loads. Zhang [4] established a translation-torsion vibration model of a compound planetary gear set, focusing on the comprehensive influence of meshing stiffness, backlash, and bearing clearance on load performance. Cooley [5] explored the modal characteristics and dynamic response prediction and then established a bending-torsion-axis-coupled model. Meanwhile, by establishing a nine-degrees-of-freedom (nine-DOF) model that considered vibration coupling between the gearbox and the shaft, Omar [6] focused on the effects of frequency, backlash, transmission error, and stiffness on nonlinear dynamic characteristics. Liu [7] analyzed the influence of gear corner contact caused by the pitch deviation in gear system dynamics. Zhang [8] introduced the influence of geometric eccentricity and shaft-bending deformation. Although many factors have been taken into account in these works, few studies have considered the nonuniformly distributed meshing force in the dynamic modeling of gear systems.

For the convenience of analysis, the meshing force is usually treated as uniformly distributed. In practical applications, however, it is nonuniformly distributed along the line of action due to machining errors, assembly errors, misalignment errors, etc. When the nonuniformly distributed meshing force is coupled with shaft deformation, the dynamic center distance, and the time-varying meshing stiffness, the transmission performance of the gear system will be seriously degraded. Therefore, the nonuniformly distributed meshing force cannot be ignored for the gear systems used in complicated working conditions, such as wind power gear systems. He [9] analyzed the effect of misalignment error on the time-varying meshing stiffness and the transmission error and found that only the errors on the plane of action significantly influence the meshing characteristics. Yang [10] studied the influence of misalignment error on the contact pattern of the gear pair and illustrated the possible harm of the nonuniformly distributed meshing force through the analysis of the contact area. Besides, Shi [11] calculated the bending moment and meshing parameter of an enhanced hypoid geared rotor system under misalignment. Zhou [12] analyzed the influence of planetary gear position variation and bearing stiffness nonlinearity on the misalignment model of planetary gears. Furthermore, Zhang [13] calculated meshing force variation based on a comprehensive gear model considering various dynamic factors and analyzed system vibration under different loads and transmission errors. Nevertheless, most of them only focused on meshing force distribution and gear wear and ignored the coupling effect of the nonuniformly distributed meshing force, the shaft deformation, the dynamic center distance, and the time-varying meshing stiffness. Besides, little attention has been paid to the vibration and stability analysis of a gear system subjected to the nonuniformly distributed meshing force.

In order to decrease gear vibration and obtain better transmission performance, many researchers have studied the nonlinear dynamic characteristics of the gear system [14,15]. Lin [16] and Xu [17] analyzed the influence of several parameters on dynamic transmission errors and provided error excitation data to control gear vibration and noise. Brancati [18] studied the nonlinear frequency-response characteristics of the spur gear system and analyzed the nonlinear vibration characteristics under multiharmonic excitation. Zhang [19] investigated the influences of gear wear and bearing clearance on the dynamic responses of a gear system, and obtained a trend of factors such as translation displacement, meshing stiffness, and bearing clearance under gradual wear. In addition, Wang [20] analyzed the nonlinear dynamic characteristics of a gear system from the perspective of system stability through bifurcation diagrams, phase diagrams, Poincaré maps, and time-domain diagrams. Xiang et al. [21,22] revealed the bifurcation and chaos of a multistage gear system. Gu [23] and Geng [24] analyzed the effect of tooth surface wear on the steady-state response of the gear system. Yang [25] studied the pressure angle and backlash changes caused by bearing deformation in order to determine whether the gear system falls into a chaotic and unstable state. Huang [26] studied the influence of excitation frequency, backlash, meshing damping ratio, and other factors on the nonlinear dynamic characteristics of a system with the aid of the stability criterion. Nonetheless, most of them ignored the nonuniformly distributed meshing force and its coupling with shaft deformation, dynamic center distance, and time-varying meshing stiffness. Therefore, it is difficult for them to provide effective guidance for designing the gear systems used in complicated working environments.

In this paper, the dynamic modeling and nonlinear analysis of a spur gear system, considering a nonuniformly distributed meshing force, are researched. This paper consists of four sections. Following the introduction, a nonuniformly distributed meshing force model is developed in Section 2 via an 18-degrees-of-freedom bending-torsion-swing-coupled dynamic model of a pair of involute spur gears. Section 3 analyzes the influence of meshing frequency, stiffness excitation, damping, and error excitation on the nonlinear dynamic characteristics of the gear system. Reasonable parameter ranges that guarantee the gear system is in a stable motion are extracted. In the last section, the conclusion is summarized.

2. Nonlinear Dynamic Modeling of Spur Gear Systems

An 18-DOF nonlinear dynamic model of a pair of involute spur gears was developed, considering the coupling of the nonuniformly distributed meshing force, the shaft bending deformation, the dynamic center distance, and the time-varying meshing stiffness. Due to the strong coupling of these factors, the gear system exhibits various complicated nonlinear dynamic characteristics.

2.1. Nonuniformly Distributed Meshing Force

The meshing force is usually treated as an ideal uniform load to facilitate an analysis of the gear system. In practical applications, however, the external load of the gear system is more or less unbalanced due to machining errors, assembly errors, misalignment errors, etc., especially for the wind turbine gear system used in complicated working conditions, which will induce the swing effect to the gears and cause complicated nonlinear dynamic responses. In this subsection, the swing effect caused by the deflection torque around the X and Y axes is taken into consideration, which can better reflect the actual dynamic behavior of a gear system.

As shown in Figure 1, plane P1P2-Q1Q2 denotes the driving and driven gear’s plane of action. According to gear meshing theory, the meshing force between the gears is in the plane of action and acts on the straight line (P1P2) along the face width. However, due to misalignment errors and the elastic deformation of low-rigidity components, the actual meshing force is no longer uniformly distributed along the face width, as shown in Figure 2.

According to the derivation of misalignment errors, the load of the unit face width can be approximately expressed, as it is in [9].

$$F(z)=F_{\max }\left[1-{\left(\frac{z}{B^{\prime }}\right)}^3\right]$$

(1)

where $F_{\max }$ denotes the maximum load, and $B^{\prime }$ is the length of the actual contact line.

Then, the meshing force $F$ can be defined as follows [9]:

$$F={\displaystyle {\int }_0^{B^{\prime }}F(z)dz=\frac{3}{4}}{B}^{\prime }{F}_{\max }$$

(2)

The maximum load can be expressed as [9]:

$$F_{\max }=C_{\gamma }{e}_m\frac{B^{\prime }}{B}$$

(3)

where $C_{\gamma }$ denotes the spring constant, and $e_m$ is the error on the plane of action.

Then the length of the contact line can be deduced according to Equations (1)–(3), i.e.,

$$B^{\prime }=\sqrt{\frac{4FB}{3C_{\gamma }{e}_m}}$$

(4)

As shown in Figure 3, the meshing force can be equivalent to the resultant force F and the deflection distance d. Thus, the deflection torque T perpendicular to the plane P₁P₂-Q₁Q₂ is obtained as:

$$T=Fd$$

(5)

where the meshing force F is derived from the Hertz contact theory: $d=\frac{B}{2}-\frac{{\displaystyle {\int }_0^{B^{\prime }}zF(z)dz}}{{\displaystyle {\int }_0^{B^{\prime }}F(z)dz}}$.

Based on the above analysis, a single-stage involute spur gear system, considering the elastic deformation of the supporting bearings and driveshafts, was researched, as shown in Figure 4. Due to the existence of the deflection moment, T, the bearing deformations are different. Therefore, it is necessary to consider two translational DOF along the X-axis and Y-axis for each supporting bearing. Meanwhile, the translational DOF along the X-axis and Y-axis and three rotational DOF around the centroid coordinate system are considered for each gear. Therefore, an 18-DOF bending-torsion-swing-coupled dynamic model is proposed in this study to depict the actual dynamic behavior of the gear system subjected to a nonuniformly distributed meshing force.

The inertial coordinate systems of the gears and the supporting bearings are ${\mathrm{O}}_i{\mathrm{X}}_i{\mathrm{Y}}_i{\mathrm{Z}}_i\left(i=1,2\right)$ and ${\mathrm{B}}_i{\mathrm{X}}_{\mathrm{b}i}{\mathrm{Y}}_{\mathrm{b}i}{\mathrm{Z}}_{\mathrm{b}i}\left(i=1,2,3,4\right)$, respectively. The global coordinate system $\mathrm{OXYZ}$ coincides with ${\mathrm{O}}_1{\mathrm{X}}_1{\mathrm{Y}}_1{\mathrm{Z}}_1$. The origin of each coordinate system is at the centroid of the gears. Besides, the axis ${\mathrm{X}}_1$ coincides with the axis ${\mathrm{X}}_2$; the axis ${\mathrm{Z}}_{bi}\left(i=1,2\right)$ coincides with the axis ${\mathrm{Z}}_1$; the axis ${\mathrm{Z}}_{bi}\left(i=3,4\right)$ coincides with the axis ${\mathrm{Z}}_2$. Meanwhile, $x_i,y_i\left(i=1,2\right)$ denotes the displacements of the driving and driven gears with respect to the origin of their own coordinate system, respectively. $x_{bi},y_{bi}\left(i=1,2,3,4\right)$ represent the displacements of the centroid of the supporting bearings.

2.2. Bending Deformation of Driveshafts

The bending deformation of the driveshaft plays an important role when subjected to a large load. In view of the swing effect caused by the nonuniformly distributed meshing force, two translational DOF are taken into account for each bearing to establish a relatively accurate dynamic model.

The schematic diagram of the driving gear’s shaft that is projected on the plane O₁YZ is shown in Figure 5. When the gear system is in a static state, $l_{bi}(i=1,2,3,4)$ denotes the distance from the centroid of the driving or driven gear to the centroid of the corresponding supporting bearing, and $L_j(j=1,2)$ denotes the length of the driving shaft or the driven shaft. Besides, along the Y-axis, the displacement of the gear centroid is $y_1$, and the ideal displacement is $y_{c1}$ when the drive shaft is regarded as a rigid body.

When the shaft is regarded as a rigid body, the gear centroid is at point A. However, the bending deformation of the driveshaft will cause the gear centroid to move to point B. Thus, one obtains the following:

$$y_{c1}=y_{b1}+\left(y_{b2}-y_{b1}\right)l_{b1}/L$$

(6)

$${\Delta }_{1}y=y_1-y_{c1}=y_1-y_{b2}{l}_{b1}/L-y_{b1}\left(L-l_{b1}\right)/L$$

(7)

where ${\Delta }_{1}y$ is the elastic deformation of the driving gear shaft along the Y-axis.

The elastic deformation of the driveshaft along the X-axis can be obtained in the same way. Thus, the elastic deformation of the driving and driven shafts at the gear centroid can be obtained as

$$\{\begin{array}{l}{\Delta }_{1}x=x_1-{\delta }_1{x}_{b2}-{\delta }_2{x}_{b1}\\ {\Delta }_{1}y=y_1-{\delta }_1{y}_{b2}-{\delta }_2{y}_{b1}\\ {\Delta }_{2}x=x_2-{\delta }_3{x}_{b4}-{\delta }_4{x}_{b3}\\ {\Delta }_{2}y=y_2-{\delta }_3{y}_{b4}-{\delta }_4{y}_{b3}\end{array}$$

(8)

where ${\delta }_i=l_{bi}/L_j$ (when $i=1,2$, $j=1$; when $i=3,4$, $j=2$).

2.3. Dynamic Center Distance

When subjected to a large load, the elastic deformation of the driveshafts and supporting bearings can cause time-varying center distances and result in the difference between the actual and theoretical line of action (LOA). As a result, the direction of the meshing force between the gears is no longer along the theoretical LOA.

As shown in Figure 6, the dashed line indicates the initial positions of the driving gear and the driven gear, while the solid line indicates their positions at an arbitrary time. The straight line PQ represents the actual dynamic LOA at an arbitrary time. L and L₀ denote the actual dynamic center distance and the theoretical center distance, respectively. Meanwhile, $\alpha$ indicates the dynamic pressure angle in the working process, and $\beta$ indicates the angle between the dynamic LOA and the X-axis.

Let the position vectors of the two gears’ centroids be R₁ and R₂; then, the center distance L of the gear pair at an arbitrary given time is

$$L=\left|{\mathbf{R}}_2-{\mathbf{R}}_1\right|=\sqrt{{\left({\mathrm{L}}_0+\Delta x\right)}^2+\Delta y^2}$$

(9)

where ${\mathbf{R}}_{\mathbf{1}}=x_1\mathbf{i}+y_1\mathbf{j}$, ${\mathbf{R}}_{\mathbf{2}}=\left(x_2+{\mathrm{L}}_0\right)\mathbf{i}+y_2\mathbf{j}$, $\Delta x=x_2-x_1$, $\Delta y=y_2-y_1$. $\mathbf{i}$ and $\mathbf{j}$ are unit vectors for the X and Y axes, respectively.

The pressure angle of the driving and driven gears can be obtained as

$$\alpha =\arccos \left(\frac{r_{b1}+r_{b2}}{L}\right)$$

(10)

where $r_{b1}$ and $r_{b2}$ are the base radii of the driving gear and the driven gear, respectively.

According to the geometric relationship of the global coordinate system, the angle between the two gears’ center line and the X-axis at an arbitrary time is

$$\beta =\arctan \left(\frac{\Delta y}{L+\Delta x}\right)$$

(11)

Therefore, taking into account the influence of elastic deformation and gear transmission errors on the gear center distance, the meshing deformation between two gears is obtained as

$$\lambda =\left(x_1-x_2\right)\sin \left(\alpha -\beta \right)-\left(y_1-y_2\right)\cos \left(\alpha -\beta \right)+r_{b1}{\theta }_{z1}-r_2{\theta }_{z2}-e(t)$$

(12)

where $e(t)$ is the gear transmission error, which causes the actual gear meshing position to deviate from the theoretical position and thus influences the meshing accuracy of the gear system.

Generally, the simple harmonic function is used to represent the gear transmission error [1], i.e.,

$$e\left(t\right)=e_0+e_1\sin (\omega t+\phi )$$

(13)

where $e_0$ and $e_1$ denote the constant value and amplitude of the transmission error, respectively; $\omega$ is the meshing frequency of the gear pair; $\phi$ is the initial phase angle, generally taking the value of 0.

Actually, the swing of the gear mainly influences the position and the deformation of the driveshaft, and its influence on the meshing force is weak. The dynamic meshing force between the gears can be described by the Hertz contact model, as shown in Figure 7. In Figure 7, $k(t)$ is the time-varying stiffness introduced in Section 2.4, and $c(t)$ is the meshing damping, i.e.,

$$c(t)=2\xi \sqrt{\frac{I_1{I}_2}{r_{b2}{}^2{I}_1+r_{b1}{}^2{I}_2}k(t)}$$

(14)

where $\xi$ is the mesh damping ratio of the gear and generally takes a value between 0.03–0.17; I₁ and I₂ are the moments of inertia, and r_b1 and r_b2 are the base circle radii of the driving and driven gear, respectively.

Besides, the backlash function $f(\lambda )$ is defined as a piecewise function [6], i.e.,

$$f\left(\lambda \right)=\{\begin{array}{l}\lambda -b\lambda >b\\ 0-b\le \lambda \le b\\ \lambda +b\lambda <-b\end{array}$$

(15)

Therefore, the dynamic meshing force is finally expressed as

$$F=k(t)f(\lambda )+c(t)\dot{f}(\lambda )$$

(16)

2.4. Time-Varying Meshing Stiffness

In order to realize the continuous transmission of power and motion, the contact ratio of a spur gear pair is between 1 and 2. Therefore, when the number of meshing teeth pairs alternates during the meshing process, the meshing stiffness will periodically change.

The parabolic method is used to solve the time-varying meshing stiffness of the gear pair. When a single pair of gear teeth are engaged in the meshing process, the meshing stiffness of the pitch point can be expressed as

$$K_c=0.8\times {10}^{3}b/q$$

(17)

$$q=0.04732+0.15551/Z_{v1}+0.25791/Z_{v2}$$

(18)

where b is the face width of the gear, $Z_{v1}$ and $Z_{v2}$ are the equivalent tooth numbers of the driving and driven gears.

The meshing stiffness can be expressed as a parabolic function with respect to X. When $X\in \left[\epsilon -1,1\right]$ (ε is the contact ratio of the gear pair), a single pair of teeth are engaged. The gear meshing stiffness can be described as

$$K_1\left(X\right)=A{X}^2+BX+C$$

(19)

When $X\in [0,\epsilon -1]$, a double pair of teeth are engaged. The meshing stiffness is

$$K_2\left(X\right)=2A{X}^2+2(A+B)X+(A+B+2C)$$

(20)

where $A=-2K_0/{\epsilon }^2$, $B=2K_0/\epsilon$, $C=K_0$; $K_0$ represents the stiffness of the meshing-in point and the meshing-out point, and generally $K_C=1.5K_0$.

Within a range of engagement, let the average meshing stiffness be $\overline{K}$, i.e.,

$$\overline{K}=A{\epsilon }^3/3+B{\epsilon }^2/2+C\epsilon$$

(21)

From Equations (17)–(21), one obtains

$$\overline{K}=8\epsilon K_c/9$$

(22)

The gear parameters used in this study are listed in

Table 1. Parameters of the gear system.

Parameter	Driving Gear	Driven Gear
Teeth number z	30	36
Mass (g)	600	850
Inertia torque ($\mathrm{kg}·{\mathrm{mm}}^2$)	600	1200
Modulus (mm)	3
Pressure angle ($°$)	20
Face width (mm)	12
Average meshing stiffness k0 (N/mm)	2.39 × 105
Tooth side clearance 2b (mm)	0.1

. The time-varying meshing stiffness of the gear system can be calculated according to Equations (19) and (20), as shown in Figure 8.

In this study, the time-varying meshing stiffness is expanded into the Fourier series in order to facilitate further analysis and calculation. The formula of Fourier transformation is

$$k\left(t\right)=k_0+{\displaystyle \sum _{i=1}^i{k}_i\sin (a\omega t+{\phi }_i)}$$

(23)

where $k_0$ is the average mesh stiffness, $k_i$ is the amplitude of the i-th harmonic term, and ${\phi }_i$ is the initial phase angle of the i-th harmonic component. $i$ is the order of the harmonic term expanded by the Fourier series, and $i=8$ in this study.

2.5. Derivation of the Dynamic Model

In order to better describe the deformation of gears, driveshafts, and bearings under a nonuniformly distributed meshing force, five DOF are introduced for each gear and two translational DOF are introduced for each bearing. Therefore, this section proposes an 18-DOF bending-torsion-swing-coupled nonlinear dynamic model based on Lagrange’s equation. Let the generalized coordinates of the driving and driven gears be $z_1={\left[x_1,y_1,x_{b1},y_{b1},x_{b2},y_{b2},{\theta }_{x1},{\theta }_{y1},{\theta }_{z1}\right]}^T$ and $z_2={\left[x_2,y_2,x_{b3},y_{b3},x_{b4},y_{b4},{\theta }_{x2},{\theta }_{y2},{\theta }_{z2}\right]}^T$.

The kinetic energy of the gear system mainly consists of the translational and rotational kinetic energy of the gears and the rotational kinetic energy of the supporting bearings. Therefore, the kinetic energy function $T_1$ of the driving gear and the kinetic energy function $T_2$ of the driven gear can be expressed as

$$T_1=\frac{1}{2}{m}_1({\dot{x}}_1{}^2+{\dot{y}}_1{}^2)+\frac{1}{2}{I}_{x1}{\dot{\theta }}_{x1}{}^2+\frac{1}{2}{I}_{y1}{\dot{\theta }}_{y1}{}^2+\frac{1}{2}{I}_{z1}{\dot{\theta }}_{z1}{}^2+\frac{1}{2}{m}_{b1}({\dot{x}}_{b1}{}^2+{\dot{y}}_{b1}{}^2)+\frac{1}{2}{m}_{b2}({\dot{x}}_{b2}{}^2+{\dot{y}}_{b2}{}^2)$$

(24a)

$$T_2=\frac{1}{2}{m}_2({\dot{x}}_2{}^2+{\dot{y}}_2{}^2)+\frac{1}{2}{I}_{x2}{\dot{\theta }}_{x2}{}^2+\frac{1}{2}{I}_{y2}{\dot{\theta }}_{y2}{}^2+\frac{1}{2}{I}_{z2}{\dot{\theta }}_{z2}{}^2+\frac{1}{2}{m}_{b3}({\dot{x}}_{b3}{}^2+{\dot{y}}_{b3}{}^2)+\frac{1}{2}{m}_{b4}({\dot{x}}_{b4}{}^2+{\dot{y}}_{b4}{}^2)$$

(24b)

where $m_i\left(i=1,2\right)$ is the mass of the driving gear and the driven gear; $m_{bj}(j=1,2,3,4)$ is the equivalent mass of the rotating body of the supporting bearings; $I_{xi},I_{yi},I_{zi}\left(i=1,2\right)$ denotes the rotational inertia of the gears.

In this study, the spring-damping model is used to calculate the potential energy and dissipation energy of the elastic components. The potential energy function and the dissipation energy function of the driving and driven gear are $U_1$, $U_2$, $V_1$ and $V_2$, respectively, i.e.,

$$U_1=\frac{1}{2}{k}_{b1}{\Delta }_1{x}^2+\frac{1}{2}{k}_{b1}{\Delta }_1{y}^2+\frac{1}{2}{k}_{bxy}(x_{b1}{}^2+y_{b1}{}^2+x_{b2}{}^2+y_{b2}{}^2)$$

(25a)

$$U_2=\frac{1}{2}{k}_{b1}{\Delta }_2{x}^2+\frac{1}{2}{k}_{b1}{\Delta }_2{y}^2+\frac{1}{2}{k}_{bxy}(x_{b3}{}^2+x_{b4}{}^2+y_{b3}{}^2+x_{b4}{}^2)$$

(25b)

$$V_1=\frac{1}{2}{c}_{b1}{\Delta }_1{\dot{x}}^2+\frac{1}{2}{c}_{b1}{\Delta }_1{\dot{y}}^2+\frac{1}{2}{c}_{bxy}({\dot{x}}_{b1}{}^2+{\dot{x}}_{b2}{}^2+{\dot{y}}_{b1}{}^2+{\dot{y}}_{b2}{}^2)$$

(25c)

$$V_2=\frac{1}{2}{c}_{b1}{\Delta }_2{\dot{x}}^2+\frac{1}{2}{c}_{b1}{\Delta }_2{\dot{y}}^2+\frac{1}{2}{c}_{bxy}({\dot{x}}_{b3}{}^2+{\dot{x}}_{b4}{}^2+{\dot{y}}_{b3}{}^2+{\dot{y}}_{b4}{}^2)$$

(25d)

where $k_{b1}$ and $c_{b1}$ are the equivalent bending stiffness and damping of the two driveshafts. $k_{bxy}$ and $c_{bxy}$ are the stiffness and damping of the four supporting bearings.

Let the input and output torques of the gear system be $T_s$ and $T_{\mathrm{c}}$, respectively, and the base radii of the driving and driven gear be $r_{b1}$ and $r_{b2}$. According to Section 2.3, the generalized forces acting on the driving and driven gears can be expressed by the column vectors $Q_1$ and $Q_2$ as

$$Q_1=\left[\begin{array}{l}-F\sin \phi \\ -F\cos \phi \\ 0\\ 0\\ 0\\ 0\\ -Fd\cos \phi \\ -Fd\sin \phi \\ T_s-F{r}_{b1}\end{array}\right],Q_2=\left[\begin{array}{l}F\sin \phi \\ F\cos \phi \\ 0\\ 0\\ 0\\ 0\\ Fd\cos \phi \\ Fd\sin \phi \\ F{r}_{b2}-T_c\end{array}\right]$$

(26)

Substitute Equations (24)–(26) into Lagrange’s equation, which is

$$\frac{d}{dt}\left(\frac{\partial T}{\partial {\dot{z}}_i}\right)-\frac{\partial T}{\partial z_i}+\frac{\partial U}{\partial z_i}+\frac{\partial V}{\partial {\dot{z}}_i}=Q_i$$

(27)

The dynamic equation of the gear system is obtained as follows:

$$\{\begin{array}{l}{m}_1{\ddot{x}}_1+k_{b1}(x_1-{\delta }_1{x}_{b2}-{\delta }_2{x}_{b1})+c_{b1}({\dot{x}}_1-{\delta }_1{\dot{x}}_{b2}-{\delta }_2{\dot{x}}_{b1})=-F\sin \phi \\ m_1{\dot{y}}_1+k_{b1}(y_1-{\delta }_1{y}_{b2}-{\delta }_2{y}_{b1})+c_{b1}({\dot{y}}_1-{\delta }_1{\dot{y}}_{b2}-{\delta }_2{\dot{y}}_{b1})=-F\cos \phi \\ m_{b1}{\ddot{x}}_{b1}-k_{b1}{\delta }_2(x_1-{\delta }_1{x}_{b2}-{\delta }_2{x}_{b1})+k_{bxy}{x}_{b1}+c_{bxy}{\dot{x}}_{b1}-c_{b1}{\delta }_2({\dot{x}}_1-{\delta }_1{\dot{x}}_{b2}-{\delta }_2{\dot{x}}_{b1})=0\\ m_{b1}{\ddot{y}}_{b1}-k_{b1}{\delta }_2(y_1-{\delta }_1{y}_{b2}-{\delta }_2{y}_{b1})+k_{bxy}{y}_{b1}+c_{bxy}{\dot{y}}_{b1}-c_{b1}{\delta }_2({\dot{y}}_1-{\delta }_1{\dot{y}}_{b2}-{\delta }_2{\dot{y}}_{b1})=0\\ m_{b2}{\ddot{x}}_{b2}-k_{b1}{\delta }_1(x_1-{\delta }_1{x}_{b2}-{\delta }_2{x}_{b1})+k_{bxy}{x}_{b2}+c_{bxy}{\dot{x}}_{b2}-c_{b1}{\delta }_1({\dot{x}}_1-{\delta }_1{\dot{x}}_{b2}-{\delta }_2{\dot{x}}_{b1})=0\\ m_{b2}{\ddot{y}}_{b2}-k_{b1}{\delta }_1(y_1-{\delta }_1{y}_{b2}-{\delta }_2{y}_{b1})+k_{bxy}{y}_{b2}+c_{bxy}{\dot{y}}_{b2}-c_{b1}{\delta }_1({\dot{y}}_1-{\delta }_1{\dot{y}}_{b2}-{\delta }_2{\dot{y}}_{b1})=0\\ I_{x1}{\ddot{\theta }}_{x1}=-Fd\cos \phi \\ I_{y1}{\ddot{\theta }}_{y1}=-Fd\sin \phi \\ I_{z1}{\ddot{\theta }}_{z1}=T_s-F{r}_{b1}\\ m_2{\ddot{x}}_2+k_{b1}(x_2-{\delta }_3{x}_{b4}-{\delta }_4{x}_{b3})+c_{b1}({\dot{x}}_2-{\delta }_3{\dot{x}}_{b4}-{\delta }_4{\dot{x}}_{b3})=F\sin \phi \\ m_2{\ddot{y}}_2+k_{b1}(y_2-{\delta }_3{y}_{b4}-{\delta }_4{y}_{b3})+c_{b1}({\dot{y}}_2-{\delta }_3{\dot{y}}_{b4}-{\delta }_4{\dot{y}}_{b3})=F\cos \phi \\ m_{b3}{\ddot{x}}_{b3}-k_{b1}{\delta }_4(x_2-{\delta }_3{x}_{b4}-{\delta }_4{x}_{b3})+k_{bxy}{x}_{b3}+c_{bxy}{\dot{x}}_{b3}-c_{b1}{\delta }_4({\dot{x}}_2-{\delta }_3{\dot{x}}_{b4}-{\delta }_4{\dot{x}}_{b3})=0\\ m_{b3}{\ddot{y}}_{b3}-k_{b1}{\delta }_4(y_2-{\delta }_3{y}_{b4}-{\delta }_4{y}_{b3})+k_{bxy}{y}_{b3}+c_{bxy}{\dot{y}}_{b3}-c_{b1}{\delta }_4({\dot{y}}_2-{\delta }_3{\dot{y}}_{b4}-{\delta }_4{\dot{y}}_{b3})=0\\ m_{b4}{\ddot{x}}_{b4}-k_{b1}{\delta }_3(x_2-{\delta }_3{x}_{b4}-{\delta }_4{x}_{b3})+k_{bxy}{x}_{b4}+c_{bxy}{\dot{x}}_{b4}-c_{b1}{\delta }_3({\dot{x}}_2-{\delta }_3{\dot{x}}_{b4}-{\delta }_4{\dot{x}}_{b3})=0\\ m_{b4}{\ddot{y}}_{b4}-k_{b1}{\delta }_3(y_2-{\delta }_3{y}_{b4}-{\delta }_4{y}_{b3})+k_{bxy}{y}_{b4}+c_{bxy}{\dot{y}}_{b4}-c_{b1}{\delta }_3({\dot{y}}_2-{\delta }_3{\dot{y}}_{b4}-{\delta }_4{\dot{y}}_{b3})=0\\ I_{x2}{\ddot{\theta }}_{x2}=Fd\cos \phi \\ I_{y2}{\ddot{\theta }}_{y2}=Fd\sin \phi \\ I_{z2}{\ddot{\theta }}_{z2}=F{r}_{b2}-T_c\end{array}$$

(28)

In order to simplify the 9th and the 18th expression in Equation (28), a new variable, $x=r_{b1}{\theta }_{z1}-r_{b2}{\theta }_{z2}-e(t)$, is introduced to represent the relative displacement of the gears along the LOA, caused by torsional vibration and gear error. Then, let $m_3={I_{z1}/r_{b1}}^2$, $m_4={I_{z2}/r_{b2}}^2$, $F_1=T_s/r_{b1}$, $F_2=T_c/r_{b2}$, and $m=m_3{m}_4/(m_3+m_4)$. The simplified equations can be obtained as, i.e.,

$$\lambda =(x_1-x_2)\sin \phi -(y_1-y_2)\cos \phi +x$$

(29)

$$m\ddot{x}+k(t)f(\lambda )+c\dot{\lambda }=\frac{F_{1}m}{m_3}+\frac{F_{2}m}{m_4}-m\ddot{e}(t)=F(t)$$

(30)

where $F(t)$ denotes the sum of the equivalent external error excitation and equivalent internal error excitation. Let $F_m=\frac{F_{1}m}{m_3}+\frac{F_{2}m}{m_4}$, and $F_h\left(t\right)=-m\ddot{e}\left(t\right)$. $F_m$ denotes the equivalent external excitation, and $F_h\left(t\right)$ denotes equivalent internal error excitation.

Meanwhile, in order to facilitate nonlinear analysis, the above dynamic equations of the gear system should be nondimensionalized. For this reason, let

$${\omega }_n=\sqrt{k_0/m}$$

(31)

where ${\omega }_n$ denotes the natural frequency of the gear system.

Then, let $\tau ={\omega }_{n}t$, $b_c$ denote the nominal scale of displacement; $X=x/b_c$, $X_i=x_i/b_c\left(i=1,2\right)$, $Y_i=y_i/b_c(i=1,2)$, $X_{bi}=x_{bi}/b_c\left(i=1,2,3,4\right)$, $Y_{bi}=y_{bi}/b_c\left(i=1,2,3,4\right)$, one obtains

$$\dot{x}=\frac{dx}{dt}=\frac{d\left(X{b}_c\right)}{d\tau }\frac{d\tau }{dt}=b_c{\omega }_n\dot{X}$$

(32)

$$\ddot{x}=\frac{d^{2}x}{d{t}^2}=\frac{d\left(b_c{\omega }_n\dot{X}\right)}{d\tau }\frac{d\tau }{dt}=b_c{\omega }_n{}^2\ddot{X}$$

(33)

$$\mathrm{ƛ}=(X_1-X_2)\sin \phi -(Y_1-Y_2)\cos \phi +X$$

(34)

Substituting Equations (32) and (33) into Equation (30), one obtains

$$\ddot{X}+\frac{c}{m{\omega }_n}\dot{\mathrm{ƛ}}+\frac{k\left(t\right)}{m{\omega }_n{}^2}f\left(\mathrm{ƛ}\right)=\frac{F\left(t\right)}{m{b}_c{\omega }_n{}^2}$$

(35)

Let $\xi =c/m{\omega }_n$, $K(t)={k(t)/m{\omega }_n}^2$, ${\omega }_h={\omega }_m/{\omega }_n$, $B=b/b_c$, $F_H\left(t\right)=F_h\left(t\right)/b_c$, and $F_M=F_m/\left(m{b}_c{\omega }_n{}^2\right)$. Then substitute them into Equation (35) to obtain

$$\ddot{X}+K(\tau )f(\mathrm{ƛ})+\xi \dot{\mathrm{ƛ}}=F_M+F_H(\tau )$$

(36)

Let $M_i=m/m_i\left(i=1,2\right)$, $K_b=k_b/k_0$, $C_b=c_b/c$, $K_{bxy}=k_{bxy}/k_0$, $C_{bxy}=c_{bxy}/c$, $K_1=k_1/k_0$, $F_e=e_1/b_c$, $M_{bi}=m/m_{bi}\left(i=1,2,3,4\right)$, and $\mathsf{\Lambda }=d{b}_c/r_1^2$. Finally, the simplified dimensionless dynamic equations of the gear system are obtained as:

$$\{\begin{array}{l}\ddot{X}+K(\tau )f(\mathrm{ƛ})+\xi \dot{\mathrm{ƛ}}=F_M+F_H(\tau )\\ {\ddot{X}}_1+M_1{K}_b(X_1-0.5X_{b2}-0.5X_{b1})+M_1{C}_b\xi ({\dot{X}}_1-0.5{\dot{X}}_{b2}-0.5{\dot{X}}_{b1})+M_{1}K(\tau )f(\mathrm{ƛ})\sin \phi +M_1\xi \dot{\mathrm{ƛ}}\sin \phi =0\\ {\ddot{Y}}_1+M_1{K}_b(Y_1-0.5Y_{b2}-0.5Y_{b1})+M_1{C}_b\xi ({\dot{Y}}_1-0.5{\dot{Y}}_{b2}-0.5{\dot{Y}}_{b1})+M_{1}K(\tau )f(\mathrm{ƛ})\cos \phi +M_1\xi \dot{\mathrm{ƛ}}\cos \phi =0\\ {\ddot{X}}_{b1}-0.5K_b{M}_{b1}(X_1-0.5X_{b2}-0.5X_{b1})+K_{bxy}{M}_{b1}{X}_{b1}-0.5C_b{M}_{b1}\xi ({\dot{X}}_1-0.5{\dot{X}}_{b2}-0.5{\dot{X}}_{b1})+C_{bxy}{M}_{b1}\xi {\dot{X}}_{b1}=0\\ {\ddot{Y}}_{b1}-0.5K_b{M}_{b1}(Y_1-0.5Y_{b2}-0.5Y_{b1})-0.5C_b{M}_{b1}\xi ({\dot{Y}}_1-0.5{\dot{Y}}_{b2}-0.5{\dot{Y}}_{b1})+K_{bxy}{M}_{b1}{Y}_{b1}+C_{bxy}{M}_{b1}\xi {\dot{Y}}_{b1}=0\\ {\ddot{X}}_{b2}-0.5K_b{M}_{b2}(X_1-0.5X_{b2}-0.5X_{b1})+K_{bxy}{M}_{b2}{X}_{b2}-0.5C_b{M}_{b2}\xi ({\dot{X}}_1-0.5{\dot{X}}_{b2}-0.5{\dot{X}}_{b1})+C_{bxy}{M}_{b2}\xi {\dot{X}}_{b2}=0\\ {\ddot{Y}}_{b2}-0.5K_b{M}_{b2}(Y_1-0.5Y_{b2}-0.5Y_{b1})-0.5C_b{M}_{b2}\xi ({\dot{Y}}_1-0.5{\dot{Y}}_{b2}-0.5{\dot{Y}}_{b1})+K_{bxy}{M}_{b2}{Y}_{b2}+C_{bxy}{M}_{b2}\xi {\dot{Y}}_{b2}=0\\ {\ddot{\theta }}_{x1}+4M_{1}K(\tau )\mathrm{ƛ}\mathsf{\Lambda }\cos \phi +4M_1\xi \dot{\mathrm{ƛ}}\mathsf{\Lambda }\cos \phi =0\\ {\ddot{\theta }}_{y1}+4M_{1}K(\tau )\mathrm{ƛ}\mathsf{\Lambda }\sin \phi +4M_1\xi \dot{\mathrm{ƛ}}\mathsf{\Lambda }\sin \phi =0\\ {\ddot{X}}_2+M_2{K}_b(X_2-0.5X_{b4}-0.5X_{b3})+M_2{C}_b\xi ({\dot{X}}_2-0.5{\dot{X}}_{b4}-0.5{\dot{X}}_{b3})-M_{2}K(\tau )f(\mathrm{ƛ})\sin \phi -M_2\xi \dot{\mathrm{ƛ}}\sin \phi =0\\ {\ddot{Y}}_2+M_2{K}_b(Y_2-0.5Y_{b4}-0.5Y_{b3})+M_2{C}_b\xi ({\dot{Y}}_2-0.5{\dot{Y}}_{b4}-0.5{\dot{Y}}_{b3})-M_{2}K(\tau )f(\mathrm{ƛ})\cos \phi -M_2\xi \dot{\mathrm{ƛ}}\cos \phi =0\\ {\ddot{X}}_{b3}-0.5K_b{M}_{b3}(X_2-0.5X_{b4}-0.5X_{b3})+K_{bxy}{M}_{b3}{X}_{b3}-0.5C_b{M}_{b3}\xi ({\dot{X}}_2-0.5{\dot{X}}_{b4}-0.5{\dot{X}}_{b3})+C_{bxy}{M}_{b3}\xi {\dot{X}}_{b3}=0\\ {\ddot{Y}}_{b3}-0.5K_b{M}_{b3}(Y_2-0.5Y_{b4}-0.5Y_{b3})-0.5C_b{M}_{b3}\xi ({\dot{Y}}_2-0.5{\dot{Y}}_{b4}-0.5{\dot{Y}}_{b3})+K_{bxy}{M}_{b3}{Y}_{b3}+C_{bxy}{M}_{b3}\xi {\dot{Y}}_{b3}=0\\ {\ddot{X}}_{b4}-0.5K_b{M}_{b4}(X_2-0.5X_{b4}-0.5X_{b3})+K_{bxy}{M}_{b4}{X}_{b4}-0.5C_b{M}_{b4}\xi ({\dot{X}}_2-0.5{\dot{X}}_{b4}-0.5{\dot{X}}_{b3})+C_{bxy}{M}_{b4}\xi {\dot{X}}_{b4}=0\\ {\ddot{Y}}_{b4}-0.5K_b{M}_{b4}(Y_2-0.5Y_{b4}-0.5Y_{b3})-0.5C_b{M}_{b4}\xi ({\dot{Y}}_2-0.5{\dot{Y}}_{b4}-0.5{\dot{Y}}_{b3})+K_{bxy}{M}_{b4}{Y}_{b4}+C_{bxy}{M}_{b4}\xi {\dot{Y}}_{b4}=0\\ {\ddot{\theta }}_{x2}-4M_{2}K(\tau )\mathrm{ƛ}\mathsf{\Lambda }\cos \phi -4M_2\xi \dot{\mathrm{ƛ}}\mathsf{\Lambda }\cos \phi =0\\ {\ddot{\theta }}_{y2}-4M_{2}K(\tau )\mathrm{ƛ}\mathsf{\Lambda }\sin \phi -4M_2\xi \dot{\mathrm{ƛ}}\mathsf{\Lambda }\sin \phi =0\end{array}$$

(37)

The 18-DOF bending-torsion-swing coupled nonlinear dynamic model proposed in this study can provide a relatively accurate description of the swing effect caused by the nonuniformly distributed meshing force, as well as the influence of the dynamic center distance caused by the deformation of the driveshafts and the supporting bearings on the meshing force.

3. Study of the Nonlinear Dynamic Characteristics of Gear System

Based on the above dynamic model considering the nonuniformly distributed meshing force, the influence of several essential parameters, i.e., meshing frequency, stiffness excitation, damping, and error excitation, on the nonlinear dynamic characteristics of the gear system was investigated from the perspective of system stability through global bifurcation diagrams, phase diagrams, Poincaré maps, and time-domain diagrams. For convenience, in the following figures, x1 indicates the dimensionless vibration displacement along the LOA, x2 indicates the dimensionless vibration velocity, and t indicates the dimensionless time.

3.1. Influence of Meshing Frequency

As one of the key parameters in the research of gear vibration, the meshing frequency is commonly used to investigate the nonlinear dynamic characteristics of a gear system. The meshing frequency, ω, mentioned here means the ratio of the meshing frequency to the natural frequency. Let the load torque T_c = 8.1 Nm, the supporting stiffness $k_{bxy}=4\times {10}^6\mathrm{N}/\mathrm{mm}$, the gear transmission error e₀ = 0.03 mm, and the meshing damping coefficient $\xi$ = 0.12. When the meshing frequency ω takes different values, the system’s bifurcation diagram, phase diagrams, Poincaré maps, and time-domain diagrams can be obtained.

As shown in Figure 9, with an increase in meshing frequency, the system exhibits an adding periodic bifurcation phenomenon accompanied by continuously increasing chaotic band intervals. When ω < 0.5, the system is in a chaotic state. When ω is approaching 0.5, the system abruptly enters the periodic motion from the chaotic motion. As shown in Figure 10, the phase diagram is a single closed curve, and the Poincaré map is a point, which proves that the system is in a stable single-period motion. When ω = 0.85, the system enters a period-doubling motion through period-doubling bifurcation and then rapidly develops into a chaotic motion. When ω = 1.69, it can be seen in Figure 11 that the phase trajectory diagram is irregular, and the Poincaré map is characterized by discrete dots, which means the system is in a chaotic state. As shown in Figure 12, when ω = 1.93, the phase diagram shows a closed curve band with a certain thickness, and the mapping points of the Poincaré map are gathered into three places, revealing that the system is in a quasiperiodic motion state. As shown in Figure 13, when ω = 2.66, the Poincaré diagram has 4 isolated points, indicating that the system is in a quadruple-periodic motion. All these characteristics reveal the complexity of gear system motion under the nonuniformly distributed meshing force.

Furthermore, when the meshing displacement jumps, the meshing frequency is usually n/2 times the system’s natural frequency, indicating the occurrence of a resonance phenomenon. In order to verify this inference, when the load torque F_m = 27 Nm, the system’s bifurcation diagram is obtained. As shown in Figure 14, when the system’s meshing frequency is almost n/2 times the system’s natural frequency, the resonance phenomenon arises.

In conclusion, when the system is subjected to a light load, the system motion exhibits adding periodic bifurcation accompanied by chaotic band intervals. By way of choosing the appropriate rotational speed, the system can fall into a stable periodic motion, thereby improving system stability. When the system is subjected to a large load, the resonance phenomenon is especially pronounced. The maximum vibration displacement may reach dozens of millimeters and thus damage the gear system. Therefore, for the gear systems subjected to a large load, the resonance phenomenon must be evaded.

3.2. Influence of Stiffness Excitation

3.2.1. Analysis of Time-Varying Meshing Stiffness

As an important internal excitation of the gear system, the time-varying meshing stiffness has a vital influence on the dynamic characteristics. The time-varying stiffness mentioned here means the dimensionless meshing stiffness, K₁, which indicates the ratio of the amplitude of the gear meshing stiffness to the mean value of the gear meshing stiffness.

Let the load torque T_c = 5.4 Nm, the rotational speed be 1500 rpm, the gear transmission error e₀ = 0.03 mm, the damping coefficient $\xi$ = 0.12, and the supporting stiffness $k_{bxy}$ be 1.5 times the average meshing stiffness. As shown in Figure 15, Figure 16, Figure 17, Figure 18, Figure 19 and Figure 20, the system’s bifurcation diagram, phase diagrams, Poincaré maps, and time-domain diagrams are obtained when the meshing stiffness takes different values.

As shown in Figure 15, the bifurcation diagram of the displacement along the LOA is obtained with an increase in K₁. When K₁ is in (0, 0.131), the system is in a stable periodic motion state. When K₁ falls into (0.16, 0.3), the quasiperiod-doubling motion gradually turns into a stable period-doubling motion. Subsequently, the stable period-doubling motion maintains until K₁ approaches 0.3. When K₁ falls into (0.3, 0.337), the period-doubling motion becomes a period-quadrupling motion through period-doubling bifurcation. Afterward, the period-quadrupling motion becomes a period-octupling motion through period-doubling bifurcation and finally evolves into a chaotic motion. When K₁ is in (0.34, 0.45), the system motion is chaotic. On the whole, with an increase in meshing stiffness, the system gradually enters a chaotic state from a single stable period through multiple period-doubling bifurcations, and the system gradually changes from a stable state to an unstable state.

The representative meshing stiffness was selected for the analysis of the dynamic characteristics. When K₁ = 0, as shown in Figure 16, the phase trajectory is in a single closed curve, and the Poincaré map has only one isolated dot, which means the system is in a periodic motion. When K₁ = 0.16, as shown in Figure 17, the Poincaré map exhibits several dots gathering into two places, which means the system falls into a quasiperiod-doubling motion. Similarly, as shown in Figure 18 and Figure 19, when K₁ = 0.25 and K₁ = 0.33, the system is in a period-doubling motion state and a period-quadrupling motion state, respectively. When K₁ = 0.45, it can be seen in Figure 20 that the Poincaré map is characterized by discrete dots, which means the system is in a chaotic state.

In conclusion, with an increase in meshing stiffness, the system’s dynamic response gradually evolves from a stable periodic motion to a chaotic motion, thereby deteriorating transmission performance. Decreasing the amplitude of the meshing stiffness tends to obtain a stable periodic motion and good stability. As a result, with the help of a reasonable contact ratio of the gears, both the bifurcation points and the chaotic motion can be effectively evaded, thereby improving system stability.

3.2.2. Analysis of Supporting Stiffness

The elastic deformation of the supporting bearing affects the meshing force, which, in turn, affects system stability. In this subsection, the influence of the supporting stiffness on the dynamic response is analyzed. Let the load torque T_c = 5.4 Nm, the rotational speed be 1500 rpm, the gear transmission error e₀ = 0.1 mm, and the damping coefficient $\xi$ = 0.1. The bifurcation diagrams of the system, with respect to the supporting stiffness K_bxy, are shown in Figure 21.

As shown in Figure 21, with an increase in K_bxy, the system motion gradually evolves from an unstable motion to a stable motion, and the maximum value of displacement along the LOA exhibits a downward trend. Compared with Figure 15, the influence of the dimensionless supporting stiffness on the displacement along the LOA is more complicated, characterized by the adding periodic bifurcation and the chaotic band. For example, when the dimensionless supporting stiffness increases from 10.07 to 13.44, the system motion turns from a quasiperiod-tripling motion into a three-band chaotic motion and then becomes a stable period-tripling motion. As shown in Figure 22b and Figure 23b, it can be seen from the mapping aggregation state of the points that when K_bxy is 3.3 and 10.3, respectively, the motion state of the system is in the quasiperiod-doubling motion and the quasi period-tripling motion. Similarly, as shown in Figure 24 and Figure 25, when K_bxy = 11, the mapping points in the Poincaré map are irregular, indicating that the system is in a chaotic motion state, while the system falls into a stable period-tripling motion because of the three isolated points shown in the Poincaré map when K_bxy = 14. Although system stability is complicated, the amplitude of the system along the LOA decreases significantly with an increase in supporting stiffness, as shown in Figure 22c, Figure 23c, Figure 24c and Figure 25c.

In order to further explore the effect of bearing supporting stiffness, the bifurcation diagrams of the displacement of the supporting bearings in the X-axis and Y-axis were obtained, as shown in Figure 26. It can be seen that the displacement of the supporting bearings, whether in the X-axis or the Y-axis, decreases with increasing supporting stiffness. When the supporting stiffness is large enough, the displacement of the supporting bearings is close to zero. If ignoring friction, the X-axis force subjected to the gear mainly comes from the dynamic center distance. Therefore, the displacement of the supporting bearings in the X-axis is very small, as shown in Figure 26a.

In conclusion, a dimensionless supporting stiffness plays an important role in the dynamic responses of the gear system. Generally, a large supporting stiffness will improve system stability and decrease vibration response.

3.3. Influence of Damping

Selecting an appropriate damping ratio helps to improve system stability and reduce vibration amplitude. Therefore, this section analyses the influence of meshing damping and supporting damping on the stability of the gear system considering a nonuniformly distributed meshing force.

3.3.1. Analysis of Meshing Damping

Let the load torque T_c = 40 Nm, and the supporting stiffness be $4\times 10^6\mathrm{N}/\mathrm{mm}$. Figure 27 shows the bifurcation diagram of the system with respect to the meshing damping ratio, $\xi$, in the range of $\xi \in \left[0,0.2\right]$.

As shown in Figure 27, with an increase in meshing damping, the vibration displacement of the system significantly decreases, and the motion state of the system gradually evolves from an unstable motion to a stable periodic motion. When $\xi <0.015$, the system is in a chaotic motion state, and its maximum vibration displacement reaches 1.5 mm. As shown in Figure 28c, the dimensionless vibration displacement is (−5, 9). Since the backlash is (−1, 1), it can be inferred that the slipping, back-to-back meshing, and impact phenomena have occurred. With an increase in meshing damping, e.g., meshing damping falls into (0.016, 0.02), the system is still in a chaotic state, but a jump phenomenon in vibration displacement arises. When $\xi >0.03$, the system changes from the chaotic motion state to the quasi-multiperiodic motion state. When $\xi =0.053$, along with the abrupt change in the vibration displacement, the system falls into the quasiperiod-doubling motion and then remains like this within a wide range. Afterward, the system’s dynamic response gradually becomes stable with the increase in meshing damping and finally converges into a stable periodic motion, as shown in Figure 29 and Figure 30. As shown in Figure 28c, Figure 29c and Figure 30c, the displacement amplitude along the LOA shows a decreasing trend as meshing damping increases.

Therefore, an increase in the meshing damping ratio, $\xi$, will make the system gradually tend toward a stable response and reduce the maximum vibration displacement of the system.

3.3.2. Analysis of Supporting Damping

Supporting damping means the ratio of supporting damping to meshing damping. Let the load torque T_c = 10.8 Nm, the rotational speed be 960 rpm, the supporting stiffness be $4\times 10^6\mathrm{N}/\mathrm{mm}$, and the damping coefficient $\xi$ = 0.12. When supporting damping, C_bxy, takes different values, the system’s bifurcation diagram, phase diagrams, Poincaré maps, and time-domain diagrams can be obtained.

As shown in Figure 31, the influence of supporting damping on the displacement along the LOA is complicated. When supporting damping falls into (1, 1.92), the system is in the four-band chaotic motion. When supporting damping is in (1.92, 10.82), the complexity of the chaos is significantly less than before, and multiple periodic motion and quasiperiodic motion occasionally arise, as shown in Figure 32. With increasing supporting damping, the system’s dynamic response gradually becomes stable, changing from the four-band chaotic motion to the 4n periodic motion. As shown in Figure 33a,b, the system is in the stable period octuple-motion state when C_bxy = 12, indicated by eight closed curves in the phase diagram and eight isolated dots in the Poincaré map. With an increase in supporting damping, the system motion changes from the periodic octuple motion to the period-quadrupling motion through inverse period-doubling bifurcations, as shown in Figure 34.

In conclusion, increasing the meshing damping and the supporting damping can improve system stability and decrease vibration response. Therefore, it is better to take measures to increase the gear system’s damping, for example, adopting high-damping materials.

3.4. Influence of Error Excitation

The transmission error is also an important factor affecting system stability. Let the load torque F_m = 5.4 Nm, the supporting stiffness be $4\times 10^6\mathrm{N}/\mathrm{mm}$, and the damping coefficient $\xi$ = 0.12. Figure 35 shows the bifurcation diagram of the system with respect to the transmission error $e_0$ in the range of $e_0\in \left[-0.2,0.2\right]\mathrm{mm}$.

As shown in Figure 35, when the transmission error is very small, i.e., (−0.115, 0.117) mm, the system is in a stable periodic motion state. Nevertheless, it can be observed in Figure 36 and Figure 37 that the vibration displacement amplitude slightly increases. With an increase in transmission error, the system motion suddenly switches to a quasiperiod-doubling motion, and then quickly changes to the stable periodic motion. Afterward, the system falls into a chaotic motion through period-doubling bifurcation, as shown in Figure 38. As shown in Figure 39, if the transmission error is large, the system is in a stable periodic motion state, where some impacts will occur. Besides, the displacement along the LOA is in (−1, 2.5) mm, as shown in Figure 39c. Since the tooth backlash is 0.1 mm, this means that slipping and back-to-back meshing have occurred.

As shown in Figure 40, although the bifurcation diagram of a large load is similar to that of a light load, its stable range is larger than the latter’s stable range.

In conclusion, when the load torque is large, the dynamic response of the gear system is not sensitive to the transmission error. Therefore, the requirements for manufacturing costs and assembly accuracy can be reduced. However, for a gear system subjected to a light load, its dynamic response is sensitive to the transmission error. Therefore, in order to evade the bifurcation points and the chaotic motion, the transmission error should be decreased by increasing manufacturing accuracy and assembly accuracy.

4. Conclusions

This paper puts forward an 18-DOF bending-torsion-swing coupled dynamic model of a pair of involute spur gears, considering a nonuniformly distributed meshing force, shaft bending deformation, dynamic center distance, and time-varying stiffness excitation. The influence of several important parameters, i.e., meshing frequency, stiffness excitation, damping, and error excitation, on the nonlinear dynamic characteristics of the gear system was researched from the perspective of system stability through global bifurcation diagrams, phase diagrams, Poincaré maps, and time-domain diagrams.

Compared to the related works in the literature review of this study, we have made considerable progress, as shown in

Table 2. Comparison with other research work.

Related Works in the Literature Review	Our Work
Bending-torsion dynamic model	Bending-torsion-swing dynamic model
Uniformly distributed meshing force	Nonuniformly distributed meshing force
6-DOF	18-DOF
Weak coupling of uniformly distributed meshing force, shaft deformation, and dynamic center distance	Strong coupling of nonuniformly distributed meshing force, shaft deformation, and dynamic center distance

. Several important conclusions were obtained:

(1)

The proposed dynamic modeling method of the spur gear system considers a nonuniformly distributed meshing force, dynamic center distance, shaft bending deformation, and time-varying stiffness excitation. This can better reflect the actual dynamic characteristics of a gear system in complicated working situations;

(2)

When the meshing frequency is close to n/2 times the system’s natural frequency, resonance will appear. An excessive meshing stiffness amplitude will cause the system to fall into a chaotic motion state, while a large supporting stiffness will improve the system’s stability. Increasing meshing damping and supporting damping can improve system stability and decrease the vibration response. Moreover, the gear system’s dynamic response is sensitive to transmission error when subjected to a light load. In order to evade the bifurcation points and the chaotic motion, the transmission error should be decreased by increasing manufacturing accuracy and assembly accuracy;

(3)

Traditionally, a gear system is designed in terms of strength theory, characterized by material selection, force analysis, stiffness verification, and fatigue-strength design. For gear systems that are used in complicated working conditions, design should be considered not in terms of only structural strength but also in terms of motion stability. By evading complicated nonlinear dynamic behaviors, such as bifurcation and chaos, the transmission performance of a gear system can effectively be improved.