Exploration of nonlinear characteristics of gear-bearing systems with friction composite effects

Gear systems exhibit strong nonlinear and non-smooth dynamics caused by time-varying meshing stiffness (TVMS), tooth friction, and backlash. The alternating engagement of the single and double teeth and the change in the direction of the friction force at the position of the pitch circle caused discontinuity of stiffness. In this study, a TVMS calculation method that considers friction was discussed, and a mechanical model of gear transmission that considers friction and tooth-side clearance was established. The actual discontinuous stiffness was evaluated using Fourier transform, and the dynamic response of the system caused by the discontinuity was compared. Through phase and Poincaré diagrams, the influence of the discontinuity in the TVMS, considering friction, on the bifurcation characteristics and dynamic meshing force was discussed, and the disengagement characteristics of the system under various periodic-motion states were investigated. The results show that a shift in the friction vector resulted in a step phenomenon in the TVMS, which became more evident with an increasing friction coefficient. For low-speed gear systems, stiffness simulation employing Fourier transform can accurately replicate the system response; however, it can obscure the stiffness step phenomenon. The range of the four-period motion generated within the chaotic behavior of the continuous system expanded from 0.003 to 0.008. The periodic motion rapidly increased the disengagement rate of the system. The disengagement rate increased rapidly from 16.8 % in a single cycle to 58.5 % in two cycles. The results of this study provide a theoretical basis for gear design.

1 Introduction

More than 100 years of research has demonstrated that the study of gear dynamics remains a challenging and interesting field of investigation. Recent research on mechanical systems has largely been focused on models involving motion-limiting constraints. In particular, theoretical challenges in gear systems involving contact and friction remain unresolved (Karagiannis and Pfeiffer, 1991; Khulief, 2013; Luo, 2012; Marques et al., 2016; Natsiavas, 2019). In all cases, the loss of continuity or smoothness leads to special bifurcations, resulting in interesting dynamics with more complex features than those observed in smooth systems (Foale and Bishop, 1994; Leine, 2006; Nordmark et al., 2009; Nusse and Yorke, 1992). The existence of backlash transforms the gear system into a typical piecewise smooth system, implying the presence of two interfaces in the state space of the system. On both sides of each interface, the system state or Jacobian matrix exhibits discontinuous characteristics. Direct numerical integration causes integration singularity, affects integration accuracy, and changes the convergence and divergence of the system. Exploring the effects of friction and backlash on the nonlinear dynamics of a gear system is of great significance (Wang, 2002).

In recent years, nonlinear vibration theory has been extensively used to establish dynamic models from different perspectives (Jiang and Liu, 2017). Numerous complex and interesting nonlinear phenomena in gear systems have been observed and studied through numerical simulations and experiments (Pan et al., 2019; Wang et al., 2012). Vaishya and Singh (2003) and Kim et al. (2010) established a gear transmission model that considers time-varying meshing stiffness (TVMS) and sliding friction, assuming a uniform load distribution between the gear teeth. The dynamic TE of a gear system was studied using the harmonic balance method and Floquet theory. However, the TVMS was represented as rectangular waves with the load evenly distributed, resulting in significant errors compared with the actual gear model. Saxena et al. (2015) proposed a TVMS calculation method based on the potential energy method considering friction and shaft misalignment and studied the effects of gear friction and shaft misalignment with cracks on the TVMS of the system. Wang et al. (2012) and Yang et al. (2022) studied the effects of friction on the periodic and chaotic responses of gear systems through numerical calculations, reporting that friction increases the components of the system's superharmonic and subharmonic responses and accelerates the transition of the system into a chaotic state. These results show that friction is an important factor affecting the gear system.

Bifurcation theory examines the structural and stability changes in a system caused by changes in its parameters. Zhiying et al. (2004) used the pseudo-fixed-point tracking method to study the changes in the periodic solution structure of a system with changes in the damping parameters and excitation frequency. This study demonstrated that this method provides better initial iterative values for studying the stability and bifurcation of periodic motion. Luo (2005), Chen et al. (2014), Li et al. (2014), and Luo and Gegg (2007) proposed a physical and mathematical model for the collision vibration of gear pairs in the dynamic domain by further refining and improving the theories of discontinuous dynamics and local singularity. Through the application of implicit function mapping dynamics and numerical solution methods, a mathematical description of a gear transmission system with clearance was developed using implicit function mapping dynamics and numerical solution methods, offering a novel analytical method. Many studies have been conducted to better describe the TVMS in gear systems (Huang et al., 2020; Margielewicz et al., 2019; Sainsot et al., 2014; Saxena et al., 2016). However, these studies assume the continuity of the TVMS. Geng et al. (2021, 2016) analyzed the variation in friction within gear systems; established a nonlinear dynamic model considering surface friction, time-varying backlash, and meshing stiffness; and studied the bifurcation characteristics of gear systems using bifurcation diagrams and spectrographs. However, these studies did not consider the influence of friction when simultaneously calculating the TVMS and dynamics.

This study aimed to explore the effect of friction on the TVMS and the dynamic response of gear systems during the meshing process. First, a nonlinear gear transmission model that incorporates sliding friction, TVMS, and backlash was developed. This model accounts for the effect of friction in both the computation of the TVMS and the analysis of the meshing forces. Subsequently, the effect of friction on the TVMS was investigated. Bifurcation diagrams of gear systems under different external excitation frequencies were used to evaluate the changes in discontinuity characteristics at bifurcation points and in chaotic spaces, and further examination was conducted on the proportional characteristics of the detachment time of the system under different motion states. Thus, the impact of actual discontinuous TVMS (DTVMS) on the vibration and noise of gear transmission systems was revealed.

2 Equations of motion of the gear system

To more accurately analyze the chattering caused by backlash in gear transmission engineering, a nonlinear dynamic model of a multi-degree-of-freedom (DOF) gear-bearing system was established, as shown in Fig. 1. To build the model, we focused on the vibration of the gear pair along the meshing line of gears, taking into account mesh stiffness, bearing support stiffness, mesh damping, clearance, and TE. Only the transverse and torsional vibrations of the two gears in the radial direction were examined when the friction between the teeth was neglected.

Figure 1 Dynamics model of a spur gear system.

In this model, T_p and T_g are the constant external torques applied to the pinion and gear, respectively. m_i, I_i, and R_i (i = 1: pinion; 2: gear) represent the mass, moment of inertia, and reference radius of the two gears, respectively. F_bi is the bearing radial pre-tension, and k_i and c_i are the equivalent bearing stiffness and damping, respectively. k_m is the meshing stiffness of the gear pair, c_m is the meshing damping of the gear, and 2b is the gear backlash. The basic movements of the system can be described as four DOFs: two reverse movements, θ₁ and θ₂, and two translational motions, y₁ and y₂. e(t) represents the static transfer error caused by manufacturing.

2.1 Time-varying mesh stiffness considering sliding friction

Friction can change the magnitudes of the tangential and radial forces, thereby altering the bending stiffness and axial compression stiffness. The frictional force f = μF is always perpendicular to the force F. The direction of the friction shifts in the pitch circle as the meshing point moves. Consequently, an improved approach for calculating the meshing stiffness of gear systems was developed that incorporates frictional forces into the TVMS calculation.

The potential energy generated through the bending, shearing, and axial compression of gear teeth can be expressed as follows using the potential energy method theory described by Jiang and Liu (2017):

$$\begin{array}{}\text{(1)}& {\displaystyle }{U}_{\mathrm{b}}={\int }_{\mathrm{0}}^d{\displaystyle \frac{M^{\mathrm{2}}}{\mathrm{2}E{I}_x}}\mathrm{d}x,\text{(2)}& {\displaystyle }{U}_{\mathrm{s}}={\int }_{\mathrm{0}}^d{\displaystyle \frac{\mathrm{1.2}{\left(F\cos {\mathit{\alpha }}_{\mathrm{1}}\right)}^{\mathrm{2}}}{\mathrm{2}G{A}_x}}\mathrm{d}x,\text{(3)}& {\displaystyle }{U}_{\mathrm{a}}={\int }_{\mathrm{0}}^d{\displaystyle \frac{F^{\mathrm{2}}{\sin }^{\mathrm{2}}{\mathit{\alpha }}_{\mathrm{1}}}{\mathrm{2}E{I}_x}}\mathrm{d}x.\end{array}$$

The distance of the meshing point along the tooth height direction is denoted by d. The definitions of d, x, and h are as shown in Fig. 2, and the values of G, I_x, and A_x can be calculated as

$$\begin{array}{}\text{(4)}& G={\displaystyle \frac{E}{\mathrm{2}(\mathrm{1}+\mathit{\nu })}};A_x=\mathrm{2}{h}_{x}L;I_x={\displaystyle \frac{\mathrm{2}}{\mathrm{3}}}{h}_x^{\mathrm{3}}L.\end{array}$$

Figure 2 Calculation model for spur gear tooth stiffness.

By using these parameters and substituting their values into Eqs. (1)–(4), the expressions for the bending mesh stiffness k_b, shear mesh stiffness k_s, and axial compressive stiffness k_a can be expressed as

$$\begin{array}{}\text{(5)}& {\displaystyle \frac{\mathrm{1}}{k_{\mathrm{b}}}}=\left\{\begin{array}{l}{\int }_{-{\mathit{\varphi }}_{\mathrm{1}}}^{{\mathit{\varphi }}_{\mathrm{2}}}\frac{\mathrm{3}\left\{\begin{array}{c}(\mathrm{1}+\cos {\mathit{\varphi }}_{\mathrm{1}}\left(\left({\mathit{\varphi }}_{\mathrm{2}}-\mathit{\alpha }\right)\sin \mathit{\alpha }-\cos \mathit{\alpha }\right))\\ -\mathit{\mu }({\mathit{\varphi }}_{\mathrm{1}}+{\mathit{\varphi }}_{\mathrm{2}}+\sin {\mathit{\varphi }}_{\mathrm{1}}\left(\left({\mathit{\varphi }}_{\mathrm{2}}-\mathit{\alpha }\right)\sin \mathit{\alpha }-\cos \mathit{\alpha }\right))\end{array}\right\}}{\mathrm{2}E{\left[\sin \mathit{\alpha }+\left({\mathit{\varphi }}_{\mathrm{2}}-\mathit{\alpha }\right)\cos \mathit{\alpha }\right]}^{\mathrm{3}}},\\ \mathit{\alpha }>{\mathit{\alpha }}_d\\ \\ {\int }_{-{\mathit{\varphi }}_{\mathrm{1}}}^{{\mathit{\varphi }}_{\mathrm{2}}}\frac{\mathrm{3}\left\{\begin{array}{c}\left(\mathrm{1}+\cos {\mathit{\varphi }}_{\mathrm{1}}\left(\left({\mathit{\varphi }}_{\mathrm{2}}-\mathit{\alpha }\right)\sin \mathit{\alpha }-\cos \mathit{\alpha }\right)\right)\\ +\mathit{\mu }\left({\mathit{\varphi }}_{\mathrm{1}}+{\mathit{\varphi }}_{\mathrm{2}}+\sin {\mathit{\varphi }}_{\mathrm{1}}\left(\left({\mathit{\varphi }}_{\mathrm{2}}-\mathit{\alpha }\right)\sin \mathit{\alpha }-\cos \mathit{\alpha }\right)\right)\end{array}\right\}}{\mathrm{2}E{\left[\sin \mathit{\alpha }+\left({\mathit{\varphi }}_{\mathrm{2}}-\mathit{\alpha }\right)\cos \mathit{\alpha }\right]}^{\mathrm{3}}},\\ \mathit{\alpha }<{\mathit{\alpha }}_d\end{array}\right.,\end{array}$$

$$\begin{array}{}\text{(6)}& {\displaystyle \frac{\mathrm{1}}{k_{\mathrm{s}}}}=\left\{\begin{array}{l}{\int }_{-{\mathit{\varphi }}_{\mathrm{1}}}^{{\mathit{\varphi }}_{\mathrm{2}}}\frac{\mathrm{1.2}(\mathrm{1}+v){\left(\cos {\mathit{\varphi }}_{\mathrm{1}}+\mathit{\mu }\sin {\mathit{\varphi }}_{\mathrm{1}}\right)}^{\mathrm{2}}\left({\mathit{\varphi }}_{\mathrm{2}}-\mathit{\alpha }\right)\cos \mathit{\alpha }}{\mathrm{2}EL\left[\sin \mathit{\alpha }+\left({\mathit{\varphi }}_{\mathrm{2}}-\mathit{\alpha }\right)\cos \mathit{\alpha }\right]},\\ \mathit{\alpha }>{\mathit{\alpha }}_d\\ \\ {\int }_{-{\mathit{\varphi }}_{\mathrm{1}}}^{{\mathit{\varphi }}_{\mathrm{2}}}\frac{\mathrm{1.2}(\mathrm{1}+v){\left(\cos {\mathit{\varphi }}_{\mathrm{1}}-\mathit{\mu }\sin {\mathit{\varphi }}_{\mathrm{1}}\right)}^{\mathrm{2}}\left({\mathit{\varphi }}_{\mathrm{2}}-\mathit{\alpha }\right)\cos \mathit{\alpha }}{\mathrm{2}EL\left[\sin \mathit{\alpha }+\left({\mathit{\varphi }}_{\mathrm{2}}-\mathit{\alpha }\right)\cos \mathit{\alpha }\right]},\\ \mathit{\alpha }<{\mathit{\alpha }}_d\end{array}\right.,\end{array}$$

$$\begin{array}{}\text{(7)}& {\displaystyle \frac{\mathrm{1}}{k_a}}=\left\{\begin{array}{l}{\int }_{-{\mathit{\varphi }}_{\mathrm{1}}}^{{\mathit{\varphi }}_{\mathrm{2}}}\frac{\left({\mathit{\varphi }}_{\mathrm{2}}-\mathit{\alpha }\right)\cos \mathit{\alpha }{\left(\sin {\mathit{\varphi }}_{\mathrm{1}}+\mathit{\mu }\cos {\mathit{\varphi }}_{\mathrm{1}}\right)}^{\mathrm{2}}}{\mathrm{2}EL\left[\sin \mathit{\alpha }+\left({\mathit{\varphi }}_{\mathrm{2}}-\mathit{\alpha }\right)\cos \mathit{\alpha }\right]},\mathit{\alpha }>{\mathit{\alpha }}_d\\ \\ {\int }_{-{\mathit{\varphi }}_{\mathrm{1}}}^{{\mathit{\varphi }}_{\mathrm{2}}}\frac{{\left({\mathit{\varphi }}_{\mathrm{2}}-\mathit{\alpha }\right)}^{\mathrm{2}}\cos \mathit{\alpha }{\left(\sin {\mathit{\varphi }}_{\mathrm{1}}\right)}^{\mathrm{2}}}{\mathrm{2}EL\left[\sin \mathit{\alpha }+\left({\mathit{\varphi }}_{\mathrm{2}}-\mathit{\alpha }\right)\cos \mathit{\alpha }\right]},\mathit{\alpha }<{\mathit{\alpha }}_d\end{array}\right..\end{array}$$

The Hertzian and fillet foundation stiffnesses do not change because of the friction. The Hertzian stiffness k_h is expressed as

$$\begin{array}{}\text{(8)}& k_{\mathrm{h}}={\displaystyle \frac{\mathit{\pi }EL}{\mathrm{4}(\mathrm{1}-{\mathit{\nu }}^{\mathrm{2}})}}.\end{array}$$

The fillet foundation stiffness is obtained by solving the gear-body-induced tooth deflection δ_f. Sainsot et al. (2014) proposed the following analytical formula:

$$\begin{array}{}\text{(9)}& \begin{array}{rl}& {\displaystyle \frac{\mathrm{1}}{k_{\mathrm{f}}}}={\mathit{\delta }}_{\mathrm{f}}={\displaystyle \frac{\mathrm{1}}{E}}{\displaystyle \frac{F}{L}}{\cos }^{\mathrm{2}}{\mathit{\alpha }}_{\mathrm{m}}\\ & [L^{\ast }\left(\mathit{\epsilon },{\mathit{\theta }}_{\mathrm{f}}\right){\left({\displaystyle \frac{u}{S_{\mathrm{f}}}}\right)}^{\mathrm{2}}+M^{\ast }\left(\mathit{\epsilon },{\mathit{\theta }}_{\mathrm{f}}\right){\displaystyle \frac{u}{S_{\mathrm{f}}}}+P^{\ast }\left(\mathit{\epsilon },{\mathit{\theta }}_{\mathrm{f}}\right)\\ & \times \left[\mathrm{1}+Q^{\ast }\left(\mathit{\epsilon },{\mathit{\theta }}_{\mathrm{f}}\right)t{g}^{\mathrm{2}}{\mathit{\alpha }}_{\mathrm{m}}\right]],\end{array}\end{array}$$

where E is Young's modulus; ν represents Poisson's ratio; L indicates the tooth width; S_f expresses the tooth thickness on the root circle occupied by the angle 2θ_f; ε = $R_{\mathrm{f}}/r_{\mathrm{in}}$, where R_f is the radius of the root circle and r_in is the radius of the hub bore; and L^∗, M^∗, P^∗, and Q^∗ are parameters depending on ε and θ_f. Further details are discussed elsewhere (Li et al., 2014). Thus, for single-tooth contact, the total effective TVMS can be expressed as

$$\begin{array}{}\text{(10)}& \begin{array}{rl}{k}_{\mathrm{m}}^j=\mathrm{1}/({\displaystyle \frac{\mathrm{1}}{k_{\mathrm{h}}}}& +{\displaystyle \frac{\mathrm{1}}{k_{\mathrm{b}\mathrm{1}}}}+{\displaystyle \frac{\mathrm{1}}{k_{\mathrm{s}\mathrm{1}}}}+{\displaystyle \frac{\mathrm{1}}{k_{\mathrm{a}\mathrm{1}}}}+{\displaystyle \frac{\mathrm{1}}{k_{\mathrm{f}\mathrm{1}}}}+{\displaystyle \frac{\mathrm{1}}{k_{\mathrm{b}\mathrm{2}}}}+{\displaystyle \frac{\mathrm{1}}{k_{\mathrm{s}\mathrm{2}}}}\\ & +{\displaystyle \frac{\mathrm{1}}{k_{\mathrm{a}\mathrm{2}}}}+{\displaystyle \frac{\mathrm{1}}{k_{\mathrm{f}\mathrm{2}}}}),\end{array}\end{array}$$

where k_h, k_b, k_s, k_a, and k_s represent the Hertzian, bending, shear, axial compressive mesh, and fillet foundation stiffnesses, respectively. Subscripts 1 and 2 denote the driving and driven gears, respectively. Similarly, the total effective TVMS is the sum of the stiffnesses of each contact tooth pair and can be expressed as

$$\begin{array}{}\text{(11)}& k_{\mathrm{m}}^n=\sum _{j=\mathrm{1}}^n\mathrm{1}/\left({\displaystyle \frac{\mathrm{1}}{k_{\mathrm{m}}^j}}\right).\end{array}$$

Most recent studies have used the Fourier transform method to simplify the discontinuous stiffness (Huang et al., 2020). To highlight the difference in system response between DTVMS and Fourier-transformed TVMS (FTTVMS), the TVMS of the gear was subjected to an n-order Fourier transform, and the FTTVMS of the system was obtained as

$$\begin{array}{}\text{(12)}& k_{\mathrm{m}}^{\mathrm{f}}=a_{\mathrm{0}}+\sum _{i=\mathrm{1}}^n{a}_i\cos \left(w_i\mathit{\tau }\right)+b_i\sin \left(w_i\mathit{\tau }\right).\end{array}$$

2.2 Equation of motion for a gear-pairing system with backlash

Numerous expressions for gear system backlash have been derived (Margielewicz et al., 2019). To better preserve the gear system state jump caused by clearance, we used a piecewise function to describe the motion equation of the gear system. Acting on the jth meshing tooth, the meshing force includes the elastic meshing force F_ki (i = 1, 2) caused by the TVMS (dynamic load of the gear transmission) and the viscous meshing force F_di (i = 1, 2) caused by meshing damping. Let the meshing force F_pi = F_ki + F_di and the tooth surface friction force be F_fi (i = 1, 2). The meshing and tooth surface friction forces are expressed as

$$\begin{array}{}\text{(13)}& \left\{\begin{array}{l}{F}_{\mathrm{pi}}^j=k_{\mathrm{m}}^j\left(t\right)\mathit{\delta }\left(t\right)+c_{\mathrm{m}}^j\dot{\mathit{\delta }}\left(t\right)\\ F_{\mathrm{fi}}^j=\mathit{\lambda }\mathit{\mu }{F}_{\mathrm{pi}}^j\end{array}\right.,\end{array}$$

where λ is the direction function of the frictional force, which is 1 when the frictional force is in the same direction as the relative motion of the gear teeth; otherwise, it transitions to −1. The force acting on the gear system for multi-tooth meshing is shown in Fig. 3. The total dynamic mesh force and total friction force can be expressed as $F_{\mathrm{pi}}^n\left(t\right)$ = $\sum _j^n{F}_{\mathrm{pi}}^j\left(t\right)$ and $F_{\mathrm{fi}}^n$ = $\sum _j^n{F}_{\mathrm{fi}}^n$ (i = 1, 2; j = 1, 2,…n), respectively. The arm of the frictional force can be calculated as

$$\begin{array}{}\text{(14)}& \left\{\begin{array}{l}{L}_{\mathrm{p}}^j\left(t\right)=R_{\mathrm{1}}\tan \left({\mathit{\theta }}_A+\text{mod}\left({\mathit{\omega }}_{\mathrm{1}}t,{\mathit{\theta }}_B\right)\right)+(j-\mathrm{1})\mathit{\pi }m,\\ L_{\mathrm{g}}^j\left(t\right)=L_{M,N}-L_{\mathrm{p}}^j\left(t\right),j=\mathrm{1},\mathrm{2},\mathrm{\dots },n\end{array}\right..\end{array}$$

In Eq. (14), α represents the pressure angle, and R_ga is the radius of the tooth-tip circle of the driven gear.

Figure 3 Meshing process of the gear pair.

Therefore, the total friction torque of the gear system can be expressed as T_fp = $\sum _j^n{F}_{\mathrm{f}\mathrm{1}}^n{L}_{\mathrm{p}}^j\left(t\right)$ and T_fg = $\sum _j^n{F}_{\mathrm{f}\mathrm{2}}^n{L}_{\mathrm{g}}^j\left(t\right)$. The motion state of the gear system can be divided into three types because of the presence of clearances. When the tooth surface is meshed, the motion equation of the system is

$$\begin{array}{}\text{(15)}& \left\{\begin{array}{l}{m}_{\mathrm{1}}{\ddot{y}}_{\mathrm{1}}+c_{\mathrm{1}}{\dot{y}}_{\mathrm{1}}+k_{\mathrm{1}}{y}_{\mathrm{1}}+F_{\mathrm{p}\mathrm{1}}^n=F_{\mathrm{b}\mathrm{1}}\\ m_{\mathrm{2}}{\ddot{y}}_{\mathrm{2}}+c_{\mathrm{2}}{\dot{y}}_{\mathrm{2}}+k_{\mathrm{2}}{y}_{\mathrm{2}}-F_{\mathrm{p}\mathrm{2}}^n=-F_{\mathrm{b}\mathrm{2}}\\ I_{\mathrm{1}}{\ddot{\mathit{\theta }}}_{\mathrm{1}}+F_{\mathrm{p}\mathrm{1}}^n{R}_{\mathrm{1}}-T_{\mathrm{fp}}=T_{\mathrm{p}}\\ I_{\mathrm{2}}{\ddot{\mathit{\theta }}}_{\mathrm{2}}+F_{\mathrm{p}\mathrm{2}}^n{R}_{\mathrm{2}}+T_{\mathrm{fg}}=-T_{\mathrm{g}}\end{array}\right..\end{array}$$

The absolute rotation equation of the system during gear disengagement is

$$\begin{array}{}\text{(16)}& \left\{\begin{array}{l}{I}_{\mathrm{1}}{\ddot{\mathit{\theta }}}_{\mathrm{1}}=T_{\mathrm{p}}\\ I_{\mathrm{2}}{\ddot{\mathit{\theta }}}_{\mathrm{2}}=-T_{\mathrm{g}}\end{array}\right..\end{array}$$

When gear transmission operates under light-load conditions, tooth back meshing or impact phenomena may occur. The corresponding system of equations is

$$\begin{array}{}\text{(17)}& \left\{\begin{array}{l}{m}_{\mathrm{1}}{\ddot{y}}_{\mathrm{1}}+c_{\mathrm{1}}{\dot{y}}_{\mathrm{1}}+k_{\mathrm{1}}{y}_{\mathrm{1}}-F_{\mathrm{p}\mathrm{1}}^n=F_{\mathrm{b}\mathrm{1}}\\ m_{\mathrm{2}}{\ddot{y}}_{\mathrm{2}}+c_{\mathrm{2}}{\dot{y}}_{\mathrm{2}}+k_{\mathrm{2}}{y}_{\mathrm{2}}+F_{\mathrm{p}\mathrm{2}}^n=-F_{\mathrm{b}\mathrm{2}}\\ I_{\mathrm{1}}{\dot{\mathit{\theta }}}_{\mathrm{1}}-F_{\mathrm{p}\mathrm{1}}^n{R}_{\mathrm{1}}+T_{\mathrm{fp}}=T_{\mathrm{p}}\\ I_{\mathrm{2}}{\dot{\mathit{\theta }}}_{\mathrm{2}}-F_{\mathrm{p}\mathrm{2}}^n{R}_{\mathrm{2}}-T_{\mathrm{fg}}=-T_{\mathrm{g}}\end{array}\right..\end{array}$$

During the gear meshing process, relative motion occurs in the gear system owing to tooth surface contact deformation. According to the gear meshing principle, the dynamic TE of the two gears along the direction of the meshing line during rotation is expressed as

$$\begin{array}{}\text{(18)}& \mathit{\delta }\left(t\right)=r_{\mathrm{1}}{\mathit{\theta }}_{\mathrm{1}}-r_{\mathrm{2}}{\mathit{\theta }}_{\mathrm{2}}+y_{\mathrm{1}}-y_{\mathrm{2}}-e\left(t\right),\end{array}$$

where e(t) = e₀ + e_asin (ω_mt) represents the static transfer error caused by manufacturing, assembly, and other factors. e₀ is the mean value, e_a is the amplitude of fluctuation, and ω_m = ω₁z₁ = ω₂z₂ represents the meshing frequency of the gear pair. By substituting the relative displacement δ and relative velocity $\dot{\mathit{\delta }}$ into Eq. (18), the equation can be simplified as

$$\begin{array}{}\text{(19)}& m_{\mathrm{e}}\left(\ddot{\mathit{\delta }}-{\ddot{y}}_{\mathrm{1}}+{\ddot{y}}_{\mathrm{2}}\right)+c\left(t\right)\sin \left(\mathit{\delta }\right)\dot{\mathit{\delta }}+K_{\mathrm{m}}\left(t\right)f_{\mathrm{m}}\left(\mathit{\delta }\right)=F_{\mathrm{m}}-m_{\mathrm{e}}\ddot{e}\left(t\right),\end{array}$$

where $J_{\mathrm{e}}=I_{\mathrm{1}}{R}_{\mathrm{2}}^{\mathrm{2}}+I_{\mathrm{2}}{R}_{\mathrm{1}}^{\mathrm{2}}$; $m_{\mathrm{e}}=I_{\mathrm{1}}{I}_{\mathrm{2}}/J_{\mathrm{e}}$; c(t)=c_mg(t); K_m(t)=k_m(t)g(t); $g\left(t\right)=\mathrm{1}+\mathit{\mu }\left(R_{\mathrm{1}}{L}_D^j\left(t\right)I_{\mathrm{2}}+R_{\mathrm{2}}{L}_{\mathrm{g}}^j\left(t\right)I_{\mathrm{1}}\right)/J_{\mathrm{e}}$; $F_{\mathrm{m}}=(I_{\mathrm{1}}{R}_{\mathrm{1}}{T}_q+I_{\mathrm{2}}{R}_{\mathrm{2}}{T}_{\mathrm{p}})/J_{\mathrm{e}}$; and f_m(δ) is a discontinuous function describing backlash, which can be expressed as

$$\begin{array}{}\text{(20)}& f_{\mathrm{m}}\left(\mathit{\delta }\right)=\left\{\begin{array}{l}\mathit{\delta }-b,\mathit{\delta }>b\\ \mathrm{0},-b\le \mathit{\delta }\le b\\ \mathit{\delta }+b_n,\mathit{\delta }<-b\end{array}\right..\end{array}$$

Here sign(δ) is a meshing state function used to represent the direction of friction:

$$\begin{array}{}\text{(21)}& \text{sign}\left(\mathit{\delta }\right)=\left\{\begin{array}{l}\mathrm{1},\mathit{\delta }>b\\ \mathrm{0},-b\le \mathit{\delta }\le b\\ -\mathrm{1},\mathit{\delta }<-b\end{array}\right..\end{array}$$

A three-DOF nonlinear differential equation was established for the gear-bearing system, considering the effect of the elastic deformation of the support and parameters. The equation of motion for the entire gear-bearing system is

$$\begin{array}{}\text{(22)}& \left\{\begin{array}{l}{m}_{\mathrm{1}}{\ddot{y}}_{\mathrm{1}}+c_{\mathrm{1}}{\dot{y}}_{\mathrm{1}}+k_{\mathrm{1}}{y}_{\mathrm{1}}+c_{\mathrm{m}}\dot{\mathit{\delta }}+k_{\mathrm{m}}{f}_{\mathrm{m}}\left(\mathit{\delta }\right)=-F_{\mathrm{b}\mathrm{1}}\\ m_{\mathrm{2}}{\ddot{y}}_{\mathrm{2}}+c_{\mathrm{2}}{\dot{y}}_{\mathrm{2}}+k_{\mathrm{2}}{y}_{\mathrm{2}}-c_{\mathrm{m}}\dot{\mathit{\delta }}-k_{\mathrm{m}}{f}_{\mathrm{m}}\left(\mathit{\delta }\right)=F_{\mathrm{b}\mathrm{2}}\\ m_e(\ddot{\mathit{\delta }}-{\ddot{y}}_{\mathrm{1}}+{\ddot{y}}_{\mathrm{2}})+c\left(t\right)\text{sign}\left(\mathit{\delta }\right)\dot{\mathit{\delta }}+K\left(t\right)f_{\mathrm{m}}\left(\mathit{\delta }\right)\\ =F_{\mathrm{m}}-m_e\ddot{e}\left(t\right)\end{array}\right..\end{array}$$

If the feature size b_c and dimensionless time τ = ω_nt are introduced, where ω_n = $\sqrt{k_{\mathrm{0}}/m_{\mathrm{e}}}$ is the natural frequency of the gear pair and k₀ is the average mesh stiffness, then the dimensionless backlash is b_n = $b/b_c$, the dimensionless relative displacement is y = $\mathit{\delta }/b_c$, y_in = $y/b_c(i=\mathrm{1},\mathrm{2})$, F_in = $F_{\mathrm{bi}}/\left(m_i{\mathit{\omega }}_n^{\mathrm{2}}{b}_c\right)(i=\mathrm{1},\mathrm{2})$, F_em = $F_{\mathrm{m}}/\left(m_{\mathrm{e}}{\mathit{\omega }}_n^{\mathrm{2}}{b}_c\right)$, F_eh = F_eΩ²sin (Ωτ), and F_e = $e_a/b_c$, Ω = ${\mathit{\omega }}_{\mathrm{m}}/{\mathit{\omega }}_n$. Substituting these expressions into Eq. (16) yields the following dimensionless equations expressed in matrix form:

$$\begin{array}{}\text{(23)}& \begin{array}{rl}& \left[\begin{array}{ccc}\mathrm{1}& \mathrm{0}& \mathrm{0}\\ \mathrm{0}& \mathrm{1}& \mathrm{0}\\ -\mathrm{1}& \mathrm{1}& \mathrm{1}\end{array}\right]\left(\begin{array}{c}{\ddot{y}}_{\mathrm{1}n}\\ {\ddot{y}}_{\mathrm{2}n}\\ \ddot{y}\end{array}\right)+\mathrm{2}\left[\begin{array}{ccc}{\mathit{\zeta }}_{\mathrm{11}}& \mathrm{0}& {\mathit{\zeta }}_{\mathrm{13}}\\ \mathrm{0}& {\mathit{\zeta }}_{\mathrm{11}}& -{\mathit{\zeta }}_{\mathrm{23}}\\ \mathrm{0}& \mathrm{0}& {\mathit{\zeta }}_{\mathrm{33}}\end{array}\right]\left(\begin{array}{c}{\dot{y}}_{\mathrm{1}n}\\ {\dot{y}}_{\mathrm{2}n}\\ \dot{y}\end{array}\right)\\ & +\left[\begin{array}{ccc}{k}_{\mathrm{11}}& \mathrm{0}& k_{\mathrm{13}}\\ \mathrm{0}& k_{\mathrm{11}}& -k_{\mathrm{23}}\\ \mathrm{0}& \mathrm{0}& k_{\mathrm{33}}\end{array}\right]\left(\begin{array}{c}{y}_{\mathrm{1}n}\\ y_{\mathrm{2}n}\\ f_{\mathrm{e}}\left(y\right)\end{array}\right)=\left(\begin{array}{c}-F_{\mathrm{1}n}\\ F_{\mathrm{2}n}\\ F_{\mathrm{em}}+F_{\mathrm{eh}}\end{array}\right),\end{array}\end{array}$$

where ζ_ii = $c_i/\left(\mathrm{2}{m}_i{\mathit{\omega }}_n\right)$, ζ_i3 = $c_{\mathrm{m}}/\left(\mathrm{2}{m}_i{\mathit{\omega }}_n\right)(i=\mathrm{1},\mathrm{2})$, ζ₃₃ = $c\left(t\right)\text{sign}\left(y\right)/\left(\mathrm{2}{m}_e{\mathit{\omega }}_n\right)$, k_ii = $k_i/\left(m_i{\mathit{\omega }}_n^{\mathrm{2}}\right)$, k_i3 = $K\left(t\right)/\left(m_i{\mathit{\omega }}_n^{\mathrm{2}}\right)(i=\mathrm{1},\mathrm{2})$, and k₃₃ = $K\left(t\right)/\left(m_e{\mathit{\omega }}_n^{\mathrm{2}}\right)$. Similarly, the general analytical dynamics framework is used to describe the spatial configuration of the system with a finite generalized coordinate set, q = $(q^{\mathrm{1}}{q}^{\mathrm{2}},\mathrm{\dots },q^n)$. The system has three DOFs, and its state vector is q = ${\left[y_{\mathrm{1}n},y_{\mathrm{2}n},y\right]}^T$.

3 Numerical simulation and discussion

Equation (23) of the new model considers the friction factor in the calculation of TVMS and the differential equation of motion. In the numerical simulation, the influence of the friction coefficient (COF) μ and external excitation frequency Ω in the new model on the bifurcation characteristics and gear disengagement ratio of the system was evaluated primarily to verify the model and obtain new results. Sample model tests were conducted using numerical data for a single-stage gear system. The TVMS and mesh damping necessitated further calculations.

3.1 Effects of friction on gear systems

The TVMS of the spur gear pair with friction was simulated using the equations presented in Sect. 2.2. The true COF is a dynamic, nonlinear, multivariable coupled function that typically fluctuates sharply between 0.03 and 0.12. Consequently, the COF was assumed to be μ = 0.1. The main parameters used in the TVMS calculations for the pinion and gear are listed in Table 1.

Table 1 Main parameters of the gear system.

Figure 4 shows the TVMS values for a healthy gear pair under two distinct conditions: with and without friction. When friction was considered, the stiffness values of the gear mesh decreased. The stiffness changed drastically when approaching and exiting the procedure. This change in meshing stiffness was expected because as the meshing point advances near the dividing circle node, the frictional force changes direction, resulting in a significant decrease in the meshing stiffness value. The vibration effect during meshing increased as the TVMS shifted. Furthermore, as the COF increased, the stepwise variations in stiffness due to friction became more evident.

Figure 4 Numerical calculation of spur gear tooth stiffness.

The TVMS of the gear system was fitted using a high-order Fourier transform. The FTTVMS was calculated using the values of a_i and b_i listed in Table 2. The relative error in the fundamental frequency was 0.2183 %.

Table 2 Main parameters of the FTTVMS.

Figure 5 shows the system FTTVMS for a COF of μ = 0.1. Compared with DTVMS, FTTVMS reduced the original discontinuous characteristics while completely masking the step characteristics of the TVMS induced by friction.

Figure 5 Comparison of the results after conducting Fourier transformation on the discontinuous TVMS.

A new excitation source for the gear system was added, considering the impact of friction on the TVMS. Specifically, the system produced new oscillations, the amplitude of which increased with an increasing COF. The source was the change in the direction of the friction in the pitch circle. However, the response of the gear system remained dominated by the excitation of alternating single gears because of the small step value of the frictional force compared with the stiffness change in single- and double-gear meshing.

3.2 Contribution of the external excitation frequency Ω

The excitation frequency Ω, which reveals the rotational gear speed, is a fundamental parameter determining the dynamic behavior of the mechanical transmission system. Because the speed of the gears is related to the external excitation frequency, understanding the relationship between the speed and system reaction requires an evaluation of the system dynamic response as the external excitation frequency varies.

Figure 6 Bifurcation diagrams of the excitation frequency Ω.

Figure 6a shows the bifurcation diagram derived from calculations considering DTVMS. Period-doubling and inverse period-doubling bifurcations occurred in the (0.956, 0.965) interval, and period-doubling bifurcations occurred in the (0.965, 0.969) interval, causing the system to undergo a double-period motion. A period-doubling bifurcation occurred in the (0.969, 0.972) interval, and the system changed from a two-fold period motion to a four-fold period motion. A period-doubling bifurcation occurred again in the (0.972, 0.976) interval, and the system changed from a four-fold period motion to an eight-fold period motion. The system exhibited evident quasiperiodic motion in the (0.976, 0.979) interval and then entered a chaotic-motion state, exhibiting a four-fold periodic motion in the (0.986, 0.989) interval. The system underwent an evolution process of period-doubling bifurcation, quasiperiodicity, and chaos under different excitation frequencies.

The dynamic response of the system was calculated using the FTTVMS. Figure 6b shows the graph corresponding to the case of using the FTTVMS. The periodic-motion characteristics and bifurcation points of the system exhibited significant differences compared with the original data. The bifurcation point at which the system changed from two-fold periodic to four-fold periodic motion changed from 0.969 to 0.974, and the evolution time from periodic to chaotic motion decreased. The chaotic motion of the system exhibited significant quasiperiodic changes in the range of (0.979, 0.984), and the four-period motion range generated in the chaotic motion of the system increased from 0.003 to 0.008.

A comparative analysis of the bifurcation diagrams showed that the discontinuous characteristics of the TVMS were masked by the Fourier transform. Considering the time-varying and discontinuous characteristics of the TVMS, the bifurcation points of the system undergo significant changes. Although the system response calculated using the FTTVMS reflects to some extent the periodic-motion characteristics of the system, the determination of bifurcation points may result in significant errors compared with those of the actual gear system. Therefore, discontinuous TVMS has theoretical and practical significance for determining the motion state of gear systems, particularly high-speed ones.

Figure 7 Dynamic response of systems under external excitation frequency Ω = 0.992.

To better observe the evolution of the dynamic characteristics of the system with a change in pinion speed n₁, Poincaré mapping was performed at the moment of gear disengagement by considering the Poincaré section each time the trajectory left the meshing surface. Its greatest advantage lies in its ability to determine the exact number of times the trajectory leaves the meshing and contact surfaces. This mapping means that at each iteration, the pinion and gear disengage; the corresponding mathematical expression is

$$\begin{array}{}\text{(24)}& {\mathit{\sigma }}_T=\left\{\left(y,\dot{y},\mathit{\theta }\right)\in R^{\mathrm{2}}\times S,\left(t=T_n,n\in N^{+}\right)\right\}.\end{array}$$

As shown in Fig. 7, the Poincaré points were visualized in a phase diagram generated by solving for 100 steady-state periods. Figure 7a and b show the results obtained using DTVMS and FTTVMS, respectively. When the external excitation frequency was set to Ω = 0.992, the system stiffness was obtained via Fourier transformation, leading to a significant deviation in the predicted system response.

Furthermore, the dynamic meshing force of the system was calculated. As illustrated in Fig. 8, the peak value of the dynamic meshing force determined using DTVMS was higher than that obtained using FTTVMS. These results highlight the critical engineering and practical implications of accounting for discontinuous TVMS to comprehensively analyze and accurately predict the dynamic behavior of gear systems.

Figure 8 Dynamic response of systems under an external excitation frequency of Ω = 0.992.

3.3 Investigation of the disengagement ratio across several motion states

In Fig. 9, the red areas indicate the disengagement state. When the system was in the P-1 or P-2 motion, the TVMS obtained using the Fourier transform fit the calculation results using the DTVMS. However, as the speed increased (i.e., the excitation frequency increased) and the system entered a multi-period response state, the response derived using the FTTVMS fitting differed significantly. The phase diagrams in Fig. 9c and d show the response characteristics with a TVMS discontinuity, revealing that the system response period derived through Fourier fitting did not converge.

Figure 9 Phase and spectrum diagrams of systems under external excitation frequencies of (a) Ω = 0.950, (b) Ω = 0.967, (c) Ω = 0.971, and (d) Ω = 0.973.

The dynamic response of the system for the quasiperiodic- and chaotic-motion states is shown in Fig. 10, which reveals that, in contrast to the periodic motion, the system-level detachment time increased during the quasiperiodic- and chaotic-motion stages. Furthermore, its chaotic reaction became more chaotic when discontinuous TVMS was considered. Its external excitation frequency was low at low speeds, and P-1 motion was observed. The dynamic properties of a system can be accurately predicted using the Fourier fitting approach. However, its stable periodic motion became more complicated as the rotational speed increased.

Figure 10 Phase and time Poincaré cross-section diagrams of the system under external excitation frequencies of (a, b) Ω = 0.973 and (c, d) Ω = 0.981.

Although the proportion of gear separation was small, the reciprocating motion of meshing and demeshing exacerbated the system shock and vibration. To investigate the adverse effects of different dynamic responses on the disengagement characteristics of the system, the disengagement proportions of period 1, period 2, period 4, period 8, the quasiperiod, and chaos (Fig. 11) were compared and analyzed under different dynamic responses. Figures 9 and 10 show that gear disengagement occurred periodically during periodic motion, whereas gear disengagement was not periodic under quasiperiodic and chaotic motions. To obtain more accurate results, we used 100 000 time periods after removing the transient to calculate the detachment proportion of the aperiodic response (quasiperiodic or chaotic).

As shown in Fig. 11, the proportion of disengagement time during the movement of P-1 was 16.8 %, and the proportions of disengagement time during P-2 and the quasiperiodic movement were 58.5 %, 58.6 %, 58.6 %, and 58.8 %. The proportion of disengagement time after the saddle junction forks of the system during the P-1 movement increased significantly, which intensified the impact and vibration of the system. However, in the subsequent period-doubling bifurcation, the proportion of disengagement time did not increase significantly. During the chaotic response (Fig. 10), the proportion of gear disengagement time in the system was 68.1 %, showing a significant increase. The disconnection ratio of the period-doubling and chaotic responses was relatively large, and the numerical value showed no significant difference; however, the disconnection ratio of the chaotic response was significantly larger than that of the period-1 response, seriously threatening the reliability and security of the system.

Figure 11 Proportion of disengagement time.

In general, the disengagement ratio of a system is affected by variations in its dynamic responses. When gear disengagement occurred, the quasiperiodic response disengagement ratio surpassed the P-1 response, whereas the disengagement time ratios of the gears in multiple periods were not noticeably modified. Under the chaotic-motion condition, the meshing state of the gears became more complex, which rapidly increased the gear disengagement ratio. Even if the system is stable during operation, a dynamic response with multiple periods is highly dangerous and should be given more consideration.

4 Conclusions

A three-DOF dynamic model of a gear-bearing system with backlash was constructed, considering the effect of friction on the TVMS and dynamic meshing force. The TVMS fluctuation law was examined for various COFs, and the system response of FTTVMS with an external excitation frequency was investigated using Fourier transform, along with the TVMS gained with the friction and alternating single and double teeth of gears. The numerical solution was obtained using the Runge–Kutta method, and bifurcation, phase, and Poincaré diagrams were combined to investigate the impact of frequency on the gear system bifurcation features. Furthermore, the disengagement ratios of gear systems in various operational scenarios were investigated. The following conclusions were drawn.

In the TVMS, friction served as a new excitation source, and the direction of the frictional force was reversed at the pitch circle, resulting in a step decrease in stiffness. This phenomenon became more pronounced as the COF increased. The alternating meshing of single and double teeth had a greater impact on the stiffness than the switching value generated by the frictional force.

The Fourier method was used to maintain stiffness, and discontinuous stiffness was found to have a significant impact on the periodic-motion characteristics of the system. As the external excitation frequency (speed) increased, the system transitioned from period-doubling bifurcation to chaos, and the discontinuous stiffness caused the system to become chaotic earlier. In particular, during high-speed operation, the influence of friction on the meshing vibrations of gear pairs could not be ignored.

The impacts of various dynamic reactions on the detachment ratio of the gear system varied. The disengagement time ratio of gears with distinct multiple periods was not significantly different; however, when only gear separation occurred, the disengagement ratio of the quasiperiodic response was greater than that of the period-1 response. Under the chaotic-motion condition, the meshing state of the gears became more complex, which readily increased the proportion of gear disengagement.

The continuous TVMS derived from the Fourier transform method can significantly reduce the simulation calculation time. However, the results are influenced by fitting accuracy and can only represent the periodic-motion changes occurring in high-speed gear transmission systems. Consequently, the calculation method presented in this paper offers improved guidance for the design of high-speed gear transmission systems.