The gearing system operates with a concentrated mass as a whole, so the concentrated mass method can be used to model the nonlinear dynamics of the gearing system. The following assumptions need to be made when modeling: the gears have only torsional motion during operation; the frictional influence of the support bearings is not considered; and the direction of the gear meshing force is always acting on the meshing line. Then the dynamics model of a pair of involute gear pairs is represented as shown in the Figure 1.

According to the kinetic model, the set of differential equations of motion can be deduced using Newton’s law as: { T 1 = I 1 θ ¨ 1 + r b 1 c h [ r b 1 θ ˙ 1 − r b 2 θ ˙ 2 − e ( t ) ] + k ( t ) r b 1 f ¯ ( x ( t ) ) ( r b 1 − r b 2 − e ( t ) ) − T 2 = I 2 θ ¨ 2 − r b 2 c h [ r b 1 θ ˙ 1 − r b 2 θ ˙ 2 − e ( t ) ] − k ( t ) r b 2 f ¯ ( x ( t ) ) ( r b 1 − r b 2 − e ( t ) ) (1) Where k(t) is the time-varying meshing stiffness; c _h is the meshing damping coefficient; e(t) is the gear meshing error; f ¯ ( x ( t ) ) is a non-analytic function describing the tooth meshing force when the tooth side clearance is available. r _bi、T _i、I _i、θ _i (i = 1,2) and other parameters are the base radius, torque, rotational inertia and torsional angular displacement of the master and driven wheel, respectively.

When the gear meshes, its meshing stiffness is not constant but is constantly changing with time, this change is the time-varying stiffness, the excitation phenomenon caused by it is called gear stiffness excitation.

Since gearing is a periodic motion, the meshing stiffness is a periodic function with certain laws. When the same gear pair is meshed, the frequency of stiffness change is the same as the frequency of internal excitation, so the meshing stiffness can be expressed in the form of Fourier series: k ( t ) = k m + ∑ j = 1 ∞ k j ⁡ cos ( j ω n t ¯ ) (2) Where, k _m refers to the average meshing stiffness; ω _h refers to the meshing frequency of the gear; k _j refers to the amplitude of the meshing stiffness.

Now taking the first-order harmonic component of the stiffness, the time-varying meshing stiffness is: k ( t ) = k m + k 1 ⁡ cos ( ω n t ¯ ) (3)

Tooth side clearance is caused by manufacturing and assembly and is a strong nonlinear factor in dynamics studies. During the motion of the gear pair, the meshing trajectory changes continuously due to the tooth-side clearance, which affects the gear meshing accuracy and life under repeated loads. Therefore, it is necessary to consider the effect of tooth side clearance when studying the dynamic characteristics of gear systems.

The gear side clearance function is a non-linear meshing force displacement function of the gear when there is side clearance, as a segmented function. Assuming a side gap of 2b, the segmentation function f (x(t)) can be expressed as [28]: f ( x ( t ) ) = f ¯ ( x ( t ) ) k m = { x ( t ) − b x ( t ) > b 0 | x ( t ) | ≤ b x ( t ) + b x ( t ) < − b (4)

If the gap is symmetric, then f ¯ ( x ( t ) ) takes the form shown in Figure 2.

When the actual meshing position of the gear profile deviates from the ideal position, displacement excitation, also called error excitation, is formed.

The error of a gear varies periodically as it meshes, so for analysis, a Fourier series expansion of the gear meshing frequency can be performed as follows [29]: e ( t ) = e m + ∑ j = 1 ∞ [ e a j ⁡ cos ( j ω n t ) + e b j ⁡ sin ( j ω n t ) ] = e m + ∑ j = 1 ∞ e j ⁡ cos ( j ω n t + φ j ) (5) Where, e _m is the average error; e _j is the harmonic component of the error amplitude.

Since the static transmission error is small, only the first-order harmonic components are considered for ease of analysis. Taking e _m = 0 when considering only the first order, the first-order harmonic component of the meshing error is: e ( t ) = e 1 ⁡ cos ( ω n t + φ 0 ) (6)

Define the displacements of the two gears along the meshing line as x₁ and x₂, and x₁ represents the active wheel displacement and x₂ represents the driven wheel displacement. Then it is obtained that: { x 1 = r b 1 θ 1 x 2 = r b 2 θ 2 (7)

The dynamic transmission error of the gear system is x d ( t ) = x 1 − x 2 (8)

Then the difference between the dynamic transmission error and the static transmission error of the system, i.e., the combined transmission error, is: x ( t ) = x 1 − x 2 − e ( t ) (9)

Dividing equations in Eq. 1 by r _b1 and r _b2, respectively, and substituting Eq. 7 to simplify gives: { F 1 = m 1 x ¨ 1 + c h [ x ˙ 1 − x ˙ 2 − e ˙ ( t ) ] + k ( t ) f ( x 1 − x 2 − e ( t ) ) − F 2 = m 2 x ¨ 2 − c h [ x ˙ 1 − x ˙ 2 − e ˙ ( t ) ] − k ( t ) f ( x 1 − x 2 − e ( t ) ) (10) Where: m _i (i = 1,2) are the equivalent masses of the master and driven wheels, expressed as m 1 = I 1 / r b 1 2 , m 2 = I 2 / r b 2 2 ; F _i(i = 1,2)are the circumferential forces acting in the direction of the engagement line for the master and driven wheels, expressed as F 1 = T 1 / r b 1 , F 2 = T 2 / r b 2 .

In order to eliminate the effect of stiffness displacement, the transmission error of Eq. 9 is substituted into Eq. 10 to organize the following: m e x ¨ + c h x ˙ + k ( t ) f ( x ) = F ¯ m + F ¯ h (11) Where: m _e is the equivalent mass of the gear; F ¯ m is the external excitation force of the gear, caused by the fluctuation of the active torque and load torque, and the fluctuation of the external excitation force is not considered here; F ¯ h is the internal excitation force of the gear, caused by the manufacturing and installation errors. m e = I 1 I 2 I 1 r b 2 2 + I 2 r b 1 2 (12) F ¯ h = − m e d 2 e d t 2 = m e e ω n 2 ⁡ sin ( ω n t + φ n ) (13)

Substituting Eqs 3, 13 into Eq. 11 yields: m e x ¨ + c h x ˙ + ( k m + k a ⁡ cos ( ω n t ) ) f ( x ) = F ¯ m + m e e ω n 2 ⁡ sin ( ω n t + φ n ) (14)

To facilitate the computational study of the dynamical properties of gear systems, the dynamical differential equations need to be dimensionless. Assume that the inherent frequency of the gear system model is: ω 0 = k m e (15)

Defining the quantization time as t = ω n t ¯ , and introducing the displacement nominal scale b _c, other variables can be defined by t ¯ and b _c. ω ¯ n = ω n ω 0 , F m = F ¯ m b c m e ω 0 2 , F h = e b c m e ω 0 2 , ξ = c h 2 m e ω 0 , ε = k a k m (16)

From this, the dimensionless analytical model of the differential equations of the dynamics model of the gear system can be derived as follows. x ¨ ( t ) + 2 ξ x ˙ ( t ) + [ 1 + ε ⁡ cos ( ω n t ) ] f ( x ( t ) ) = F m + F h ω n 2 ⁡ sin ( ω n t + φ h ) (17)

When the gear system is under light load, the relevant parameters are taken as follows: error amplitude F _h = 0.3; average load F _m = 0.1; time-varying meshing stiffness coefficient s = 0.1; tooth side clearance b = 0.5; dimensionless frequency ω = 1; the range of meshing damping ratio is 0.01–0.2. The bifurcation diagram and the maximum Lyapunov exponent diagram of the gear system with the variation of the damping ratio are shown in Figure 3.

By observing Figure 3, it can be found that with the gradual increase of the damping ratio, the gear system first shows a complex chaotic motion until the damping ratio is 0.11, then the system enters the proposed periodic motion, and then enters a short 4 times the proposed periodic motion, and the damping ratio continues to increase, the gear system enters a stable 2 times the periodic motion. Therefore, when the mesh damping ratio gradually increases from 0.01 to 0.2, the overall gear system undergoes a change of “chaotic motion—2 times anthropomorphic periodic motion—4 times anthropomorphic periodic motion—2 times anthropomorphic periodic motion”, and its corresponding maximum Lyapunov exponent graph also undergoes an overall change of “positive—zero—negative".

When the damping factor ξ = 0.05, the dynamics of the gear system is shown in Figure 4, the time domain diagram has no obvious periodicity, the phase plane diagram does not repeat and fills a closed area, the frequency spectrum is a continuous frequency band, and the Poincare section diagram consists of patches of dense points, so the system is in a complex chaotic motion at this time. When the damping factor ξ = 0.11, the dynamics of the gear system is shown in Figure 5, and its time domain diagram has regular peaks and presents a certain periodicity, the phase diagram is presented as a closed curve band, the spectrum diagram is discrete, and the Poincare section diagram is composed of 2 point sets, so the system is in a 2-fold proposed periodic motion at this time.

When the damping factor ξ = 0.113, the dynamics of the gear system is shown in Figure 6, whose time domain diagram is a periodic curve with regular peaks, the phase diagram shows a closed curve, the spectrum diagram is discrete, and the Poincare cross-sectional diagram consists of several point sets and the number of point sets is 4. Therefore, the system is in a 4-fold proposed periodic motion at this time. When the damping factor ξ = 0.16, the motion characteristics of the gear system are shown in Figure 7. The time domain diagram has certain periodicity and the peak value remains stable, the phase diagram shows a closed curve, the spectrum diagram is discrete, and the Poincare section diagram consists of 2 points, so the system is in a 2-fold periodic motion at this time.

The kinematic characteristics of the system presented in Figures 4–7 are consistent with the overall kinematic characteristics that vary with the damping ratio shown in Figure 3.

The results show that the variation of the damping ratio coefficient has a large effect on the nonlinear dynamic characteristics of the system when the gear system is under light load. Therefore, when analyzing and verifying the gear system under light load, the stability of the system can be enhanced by increasing the damping ratio coefficient to avoid the chaotic motion.

When the gear system is under heavy load, the relevant parameters of the system are taken as follows: error amplitude F _h = 0.1; average load F _m = 0.1; time-varying meshing stiffness coefficient s = 0.1; tooth side clearance b = 0.5; dimensionless frequency ω = 1; the range of meshing damping ratio is 0.01–0.2. The bifurcation diagram and the maximum Lyapunov exponent diagram of the gear system with the variation of the damping ratio are shown in Figure 8.

By observing Figure 8, it can be found that with the gradual increase of the damping ratio, the gear system first enters a brief chaotic motion state, followed by alternating single-fold periodic motion and chaotic motion states. Until the damping ratio increases to 0.028, the system enters into a brief triple anthropomorphic periodic motion, and as the damping ratio continues to increase, the system changes to a single periodic motion state and remains stable. The maximum Lyapunov exponent undergoes a general change of “positive-zero-negative-zero-positive-zero-negative”, which is consistent with the change of motion shown in the bifurcation diagram.

When the damping ratio ξ = 0.01, the motion characteristics of the gear system are shown in Figure 9, the time domain diagram has no obvious periodicity, the peak is not regular, the phase plane diagram does not repeat and fills a closed area, the spectrum diagram is a continuous curve, and the Poincare section diagram is composed of a piece of dense points, so the system is in a complex chaotic motion state at this time. When the damping ratio coefficient increases to ξ = 0.021, the motion characteristics of the system are shown in Figure 10, the time domain diagram shows obvious periodicity and the peak has regularity, the phase diagram is a closed curve band, the spectrum diagram is a discrete curve, and the Poincare section diagram consists of 1 point, so the system is in a single periodic motion at this time.

When the damping ratio coefficient continues to increase to ξ = 0.032, the motion characteristics of the system are shown in Figure 11. The time domain diagram is a periodic curve with peak law, the phase diagram is a closed circular curve, the spectrum diagram is a curve with three peaks, and the Poincare cross-sectional diagram consists of three point sets, so the system is in a 3-fold proposed periodic motion at this time. When the damping ratio coefficient ξ = 0.16, the motion characteristics of the system are shown in Figure 12. The time domain diagram shows an obvious periodicity, and the peak is kept stable, the phase diagram is a closed curve, the spectrum diagram is a discrete curve with a single peak, and the Poincare cross-section diagram consists of a single point, so the system is in a stable single-fold periodic motion at this time.

The changes in the kinematic state of the gear system shown in Figures 9–12 are consistent with the effect of the damping ratio coefficient ξ on the system stability obtained in Figure 8, verifying the accuracy of the analysis process and results.

The results show that the variation of the damping ratio coefficient of the gear system under heavy load condition has less effect on the stability of the system.