Numerical modeling and dynamic characteristics study of coupling vibration of multistage face gearsplanetary transmission

A novel transmission with multistage face gears as the core component achieves variable speeds via differential gear shifting. Single/multistage coupled vibration models have been established in this study to derive the coupled vibration equation in order to accurately solve the load distribution between the meshing teeth and the vibration shock between the shifting stages in the transmission process, improve the transmission smoothness of the face gears during the shifting processing, suppress the resonance of face gears meshing, reduce the noise, and optimize the power transmission performance. The characterization relationships of the key parameters such as equivalent mass, rotational inertia, equivalent mesh stiffness, support stiffness, and meshing damping coefficient to dynamic characteristics were investigated. The linear and nonlinear dynamic characteristics of coupled vibration differential equations were solved. The influence rules of factors such as integrated transmission error, dynamic load, tooth surface friction, loading speed, and load on the transmission system were analyzed. The results of the study provide a theoretical basis for the expansion of field of application of transmission devices.

1 Introduction

The mechanism of multistage face gears integrates the general advantages of most planetary transmission devices, and has a novel structure and flexible arrangement, which can be widely applied to high-speed or heavy-load fields such as automobiles, ships, and airplanes. The vibration characteristics of the device have a significant impact on the transmission stability and reliability. A significant research work has been conducted by domestic and foreign scholars to deeply understand the gear dynamic performance. Sakaridis et al. (2019) investigated the effect of the individual tooth inertia on spur gear dynamics. A lumped element model was developed, in which tooth and gear inertial properties were decoupled, with the incorporation of the meshing position and load dependent meshing stiffness and backlash (Sakaridis et al., 2019). Tatar et al. (2019) developed six degrees of freedom dynamic model of a planetary geared rotor system with equally spaced planets by considering the gyroscopic effects. The dynamic model was created using a lumped parameter model of the planetary gearbox and a finite element model of the rotating shafts using Timoshenko beams (Tatar et al., 2019). Liu et al. (2019) investigated an excitation pattern of the time-varying instantaneous center of the noncircular face gear drive system, wherein a torsional-lateral-axial coupled dynamic model of the noncircular face gear system was presented under the compound parametric excitations of the time-varying instantaneous center and mesh stiffness (Liu et al., 2019). Dadon et al. (2019) focused on the fault detection and characterization based on the vibrations in a spur gear transmission. The faults were thoroughly analyzed to understand the fault manifestation in the vibration signature and find condition indicators that were robust and sensitive to the existence and severity of the fault (Dadon et al., 2019). Yu and Mechefske (2019) defined two types of helical gear pairs based on the relationship between the transverse contact ratio and overlap the contact ratio, and proposed an improved analytical model using the slicing principle for the calculation of single mesh stiffness of helical gears (Yu and Mechefske, 2019). Ren et al. (2018) presented a novel and generalized bending-torsional-axial coupled dynamic model of a herringbone planetary gear train (HPGT) to investigate the dynamic floating performances with the application of the lumped-parameter approach. The variable step Runge–Kutta algorithm was utilized to compute the dynamic responses of the HPGT system (Ren et al., 2018). Saxena et al. (2018) performed an experimental study to discover the effect of gear pair contact on the modal behavior of an actual geared rotor system mounted on the rolling element bearings. The frequency response functions of the uncoupled and coupled geared rotor systems were measured to uncover the effect of gear pair contact on the natural frequencies of the system (Saxena et al., 2018). Dong and Hu (2018) established the bend-torsion coupling dynamics load-sharing model of the helicopter face gear split torque transmission system by using the concentrated quality standards, to analyze the dynamic load-sharing characteristics, which provided a theoretical basis and data support for its dynamic performance optimization design (Dong and Hu, 2018). Atanasovska and Hedrih (2018) used a model of a central centric collision of two fictive rolling disks, with radii equal to the radii of pitch diameters of coupled gears, in order to describe the vibro-impact dynamics of spur gears. A set of transmission error equations for the calculation of the disturbance angular velocity of the pinion was developed; the transmission error had a role of excitation of vibro-impact vibrations in the gear tooth contact (Atanasovska and Hedrih, 2018). L. N. Zhang et al. (2017) proposed a dynamic model of an N-stage power-split planetary gear set. Natural frequencies and vibration modes were obtained by solving the eigenvalue problem governed by the corresponding equations (L. N. Zhang et al., 2017). Chowdhury and Yedavalli (2017) developed analytical model for the linear vibration analysis of a pair of coupled spur geared shaft system. The model was a hybrid discrete-continuous one, wherein the gears were modeled as rigid disks mounted on the elastic shafts having transverse as well as torsional flexibilities, and supported by the rigid bearings. The non-dimensional governing equations along with the natural boundary conditions were developed using the Hamilton's principle (Chowdhury and Yedavalli, 2017). A. Q. Zhang et al. (2017) established the coupled lateral-torsional-axial vibration dynamical model of the planetary gear system under a generalized coordinate system using the shafting element method. The results indicate that the flexibility of ring-gear had a great effect on the natural frequency (A. Q. Zhang et al., 2017). Hmida et al. (2017) investigated a mechanical system composed of an electric motor, an elastic coupling, and one stage spur gear to analyze the dynamic behavior of elastic coupling, its influence on the natural frequencies, the vibrations, and the response of the system. A lumped parameter model of the corresponding system was developed to identify the natural frequencies and vibration modes (Hmida et al., 2017). Wu et al. (2017) set up a nonlinear vibration model of the gear system with multiple degrees of freedom, multi-gap, variable parameter, and bending-torsion coupling by using the mass-concentrated method for NN-type small tooth-difference planetary gear system. The nonlinear vibration differential equations deduced by Lagrange equation were solved with the fourth-order Runge–Kutta-Fehlberg method. MATLAB program was adopted to systematically analyze the vibration displacement and velocity response of the gears (Wu et al., 2017). Fernandez-Del-Rincon et al. (2017) presented a model that could simultaneously consider the internal excitations due to variable meshing stiffness (including the coupling among successive tooth pairs in contact, non-linearity linked with the contacts between surfaces and dissipative effects). The model could also simulate gear dynamics in a realistic torque dependent scenario. The proposed model combined a hybrid formulation for calculation of meshing forces with a non-linear variable compliance approach for bearings. The meshing forces were obtained by means of a double approach which combined numerical and analytical aspects (Fernandez-Del-Rincon et al., 2017). Bouslema et al. (2017) aimed at the application of a substructure methodology, based on the Frequency Response Function (FRF) simulation technique, to analyze the vibration of a stage reducer connected by a rigid coupling to a planetary gear system. The computation of the vibration response was achieved by using the FRF-based substructuring method. First of all, the two subsystems were separately analyzed and their FRF were obtained. Then, the coupled model was analyzed indirectly using the substructuring technique (Bouslema et al., 2017). Zhang et al. (2016) first presented a translational-rotational dynamic model of a two-stage closed-form planetary gear set under consideration of the rotational and translational displacements to investigate the dynamic response and avoid resonance (Zhang et al., 2016). Xiao et al. (2016) proposed a coupling dynamic model for multi-stage planetary gears train (PGT) based on the gearing theory and Lagrange equation. The modified transverse-torsional model was established in the rotating Cartesian coordinates by using the lumped-parameter method; thus, the model was more accurate than the purely torsional model for describing the physical dynamics (Xiao et al., 2016). Liu et al. (2016) performed a systematic analysis of the dynamics of a helical face-gear system with 8 degrees of freedom under complex excitation. The nonlinear dynamic system was solved using the Runge–Kutta method. The bifurcation and dynamic load characteristics of the system were identified from a series of diagrams. The effect of multi-factor on the bifurcation diagrams was also analyzed (Liu et al., 2016). Cai and Lin (2017) presented a generalized nonlinear dynamic model based on Lagrange Bond graphs to study the nonlinear dynamic characteristics of a curve-face gear drive, and indicated that the vibration response of the curve-face gear is more complex than that of the face gear due to the influence of external excitation (Cai and Lin, 2017). Lin et al. (2015) established a nonlinear dynamic model of the twisting vibration of the orthogonal curve-face gear transmission, and presented the nonlinear vibration differential equation. Using the fourth-order Runge–Kutta numerical integral method to solve the equation with MATLAB, the dynamic response of orthogonal curve-face gear transmission was obtained (Lin et al., 2015). Peng and DeSmidt (2015) established a structural dynamics model for the torsional vibration of gearboxes containing a face-gear drive by considering the flexibilities of gear teeth and transmission shafts. The face-gear mesh-induced parametric instability phenomena was numerically explored via Floquet theory for various shaft characteristics and system inertia distributions (Peng and DeSmidt, 2015). Hu et al. (2015) proposed a fourteen degree-of-freedom (DOE) lumped parameter dynamic model considering the coupled translation-rotation vibration. The jump phenomenon, periodic window, doubling periodic bifurcation and chaotic behavior of the system were observed and a path to chaotic motion through doubling period bifurcation motion was also observed (Hu et al., 2015).

These results provide good references and background for the current study. The multistage face gear transmission mechanism had steady-state linear and strong nonlinear characteristics because of the effect of the factors such as structural parameters and excitation. This study aims to establish a suitable coupled vibration model that is solved via Runge–Kutta numerical method. The related dynamic performance was studied from the perspective of linear and nonlinear dynamics, and the effects of working conditions, geometry, and corresponding physical parameters on the dynamic characteristics were deeply analyzed.

2 Establishment of multistage face gears coupled vibration model

The determination and evaluation of the vibration characteristics are primary content in the study of gear dynamic characteristics. Owing to the effects of error excitations such as design, manufacturing, and assembly, the vibration is easily generated during the transmission process, which affects the transmission efficiency and performance of the system, and even leads to broken teeth in serious cases. Therefore, it is necessary to establish a single/multistage coupled vibration model of the key components of the transmission system to derive the expressions of the primary characterization parameters. Hence, the mechanism, size, and properties of the main excitation of the transmission system were analyzed, the vibration responses of the internal and external excitation were determined. The vibration influence mechanisms of the transmission system were also suggested, thus, an optimization design method of high-efficiency transmission system has been proposed.

2.1 Single/Multistage coupled vibration model

The face gear transmission system is composed of cylindrical gears and face gears meshing with each other. The meshing dynamics models under elastic support conditions are shown in Figs. 1 and 2. According to the parameters concentrated method, the cylindrical gear and face gear are regarded as the concentrated mass and concentrated inertia. The support axle is considered a massless rigid body, the bearings and the elastic supports are equivalent to the springs and the dampers. In the model, there were frictions between the stages of face gears, and dampers were added for the simulation. In the transmission process, there were normal dynamic meshing forces on the meshing teeth surface. According to the loading characteristics of face gears, the radial forces can be negligible. Similarly, the axial forces of the cylindrical gears are not working. K_cx and K_cz are the stiffness coefficients of the cylindrical gear in the x and z directions, respectively. K_fx and K_fy are the stiffness coefficients of the face gear in the x and z directions, respectively. K_h is the stiffness coefficient of the gear pair. C_cx and C_cy are the damping coefficients of the cylindrical gear in the x and z directions, respectively. C_fx and C_fy are the damping coefficients of the face gear in the x and y directions, respectively. C_h is the damping coefficient of the gear pair. e_h is the static integrated transmission error of the gear pair. b_h is the backlash of the gear pair. μ is the time-varying friction coefficient of the gear pair. θ_cy and θ_fz are the torsional angular displacements of cylindrical gear and face gear, respectively. T_c and T_f are the torques of cylindrical gear and face gear, respectively. C_f12 and C_f23 are the damping coefficients of friction and backlash between the stages of multistage face gears.

Figure 1 Kinematics model of single-stage face gear.

Figure 2 Kinematics model of multistage face gears.

In the transmission process, there was a normal dynamic meshing force between the two tooth surfaces. There were circumferential force and radial force on the cylindrical gear, and the circumferential component force and the axial force were also generated on the face gears. Therefore, the coupled vibration model contains seven degrees of freedom:

$$\begin{array}{}\text{(1)}& {\left\{x_{\mathrm{c}},z_{\mathrm{c}},{\mathit{\theta }}_{\mathrm{c}y},x_{\mathrm{f}},z_{\mathrm{f}},{\mathit{\theta }}_{\mathrm{f}x},{\mathit{\theta }}_{\mathrm{f}z}\right\}}^{\mathrm{T}}\end{array}$$

Among these:

1. Bending vibration: The translations of cylindrical gear in x axis and z axis are x_c and z_c. The translation of face gear in x axis is x_f.

2. Torsional pendulum vibration: The torsional pendulum of face gear in x axis is θ_fx.

3. Torsional vibration: The rotation of cylindrical gear in y axis is θ_cy. The rotation of face gear in z axis is θ_fz .

4. Axial vibration: The translation of face gear in z axis is z_f.

When the unilateral constraint condition with impenetrable contact was introduced, the corresponding dynamic equation is:

$$\begin{array}{}\text{(2)}& \left\{\begin{array}{l}\mathbf{M}\ddot{q}+\mathbf{C}\dot{q}+\mathbf{K}q+{\mathit{\varsigma }}_q^{\mathrm{T}}\mathit{\lambda }=\mathbf{F}+\mathbf{Q}\\ \mathit{\varsigma }(q,t)=\mathrm{0}\end{array}\right.\end{array}$$

In Eq. (2), M, C and K are the generalized mass matrix, generalized damped matrix, and generalized stiffness matrix, respectively. ς is the constraint equation, λ is the Lagrange multiplier, Q is the generalized force matrix, and F is the contact collision force.

2.2 Modeling of coupled vibration equations

According to the force analysis of the teeth, the normal dynamic meshing force F_n, and the component forces F_x and F_z along the coordinate axes x and z are:

$$\begin{array}{}\text{(3)}& \left\{\begin{array}{l}{\mathbf{F}}_{\mathrm{n}}=k_{\mathrm{m}}\left(t\right)f\left(d_{\mathrm{n}}\right)+c_{\mathrm{m}}{\dot{d}}_{\mathrm{n}}\\ {\mathbf{F}}_x={\mathbf{F}}_{\mathrm{n}}\cos {\mathit{\alpha }}_{\mathrm{n}}\\ {\mathbf{F}}_z={\mathbf{F}}_{\mathrm{n}}\sin {\mathit{\alpha }}_{\mathrm{n}}\end{array}\right.\end{array}$$

In Eq. (3), k_m(t) is the time-varying meshing stiffness, f(d_n) is the backlash function, c_m is the meshing damping, and d_n is the relative displacement in the normal direction.

The face gear is excited by vibration and error during the transmission process to generate relative displacement in the normal direction of the meshing points.

$$\begin{array}{}\text{(4)}& \begin{array}{rl}{d}_{\mathrm{n}}& =(x_{\mathrm{c}}+r_{\mathrm{cm}}{\mathit{\theta }}_{\mathrm{c}y})\cos {\mathit{\alpha }}_{\mathrm{n}}+z_{\mathrm{c}}\sin {\mathit{\alpha }}_{\mathrm{n}}-(x_{\mathrm{f}}+r_{\mathrm{fm}}{\mathit{\theta }}_{\mathrm{f}z})\\ & \cos {\mathit{\alpha }}_{\mathrm{n}}-(z_{\mathrm{f}}+r_{\mathrm{fm}}{\mathit{\theta }}_{\mathrm{f}x})\sin {\mathit{\alpha }}_{\mathrm{n}}-e_{\mathrm{n}}\left(t\right)\end{array}\end{array}$$

In Eq. (4), r_cm and r_fm are the meshing point radii of the cylindrical gear and face gear, respectively. α_n are the normal pressure angles of meshing points. e_n(t) is the transmission error in the normal direction of the face gear.

Dynamic differential equations were established for each vibration direction of face gear transmission system:

$$\begin{array}{}\text{(5)}& \left\{\begin{array}{l}{\mathbf{M}}_{\mathrm{c}}{\ddot{x}}_{\mathrm{c}}+c_{\mathrm{c}x}{\dot{x}}_{\mathrm{c}}+k_{\mathrm{c}x}{x}_{\mathrm{c}}=-{\mathbf{F}}_x-\mathit{\mu }{\overline{f}}_{\mathit{\mu }}{\mathbf{F}}_{\mathrm{n}}\sin {\mathit{\alpha }}_{\mathrm{n}}\\ {\mathbf{M}}_{\mathrm{c}}{\ddot{z}}_{\mathrm{c}}+c_{\mathrm{c}z}{\dot{z}}_{\mathrm{c}}+k_{\mathrm{c}z}{z}_{\mathrm{c}}=-{\mathbf{F}}_z-\mathit{\mu }{\overline{f}}_{\mathit{\mu }}{\mathbf{F}}_{\mathrm{n}}\cos {\mathit{\alpha }}_{\mathrm{n}}\\ {\mathbf{J}}_{\mathrm{c}y}{\ddot{\mathit{\theta }}}_{\mathrm{c}y}={\mathbf{T}}_{\mathrm{c}}-{\mathbf{F}}_{\mathrm{n}}{r}_{\mathrm{cm}}+\mathit{\mu }{\overline{f}}_{\mathit{\mu }}{L}_A{\mathbf{F}}_{\mathrm{n}}\sin {\mathit{\alpha }}_{\mathrm{n}}\\ {\mathbf{M}}_{\mathrm{f}}{\ddot{x}}_{\mathrm{f}}+c_{\mathrm{f}x}{\dot{x}}_{\mathrm{f}}+k_{\mathrm{f}x}{x}_{\mathrm{f}}={\mathbf{F}}_x+\mathit{\mu }{\overline{f}}_{\mathit{\mu }}{\mathbf{F}}_{\mathrm{n}}\sin {\mathit{\alpha }}_{\mathrm{n}}\\ {\mathbf{M}}_{\mathrm{f}}{\ddot{z}}_{\mathrm{f}}+c_{\mathrm{f}z}{\dot{z}}_{\mathrm{f}}+k_{\mathrm{f}z}{z}_{\mathrm{f}}={\mathbf{F}}_z+\mathit{\mu }{\overline{f}}_{\mathit{\mu }}{\mathbf{F}}_{\mathrm{n}}\cos {\mathit{\alpha }}_{\mathrm{n}}\\ {\mathbf{J}}_{\mathrm{f}x}{\ddot{\mathit{\theta }}}_{\mathrm{f}x}+c_{{\mathit{\theta }}_{\mathrm{f}x}}{\dot{\mathit{\theta }}}_{\mathrm{f}x}+k_{{\mathit{\theta }}_{\mathrm{f}x}}{\mathit{\theta }}_{\mathrm{f}x}={\mathbf{F}}_z{r}_{\mathrm{fm}}-\mathit{\mu }{\overline{f}}_{\mathit{\mu }}{\mathbf{F}}_{\mathrm{n}}\sin {\mathit{\alpha }}_{\mathrm{n}}\\ {\mathbf{J}}_{\mathrm{f}z}{\ddot{\mathit{\theta }}}_{\mathrm{f}z}={\mathbf{F}}_{\mathrm{n}}{r}_{\mathrm{fm}}-{\mathbf{T}}_{\mathrm{f}}+\mathit{\mu }{\overline{f}}_{\mathit{\mu }}{r}_{\mathrm{fm}}{\mathbf{F}}_{\mathrm{n}}\cos {\mathit{\alpha }}_{\mathrm{n}}\end{array}\right.\end{array}$$

In Eq. (5), M_c and M_f are the concentrated mass of cylindrical gear and face gear, respectively. k_ij and $c_{ij}(i=c,f;j=x,z)$ are the support stiffness and damping of the cylindrical gear or face gear along the x and z axes, respectively. J_cy is the rotational inertia of the cylindrical gear along the y axes. J_fx and J_fz are the rotational inertias of face gear in the x and z axes, respectively. T_c is the input torque acting on the cylindrical gear. T_f is the load torque acting on face gear. $c_{{\mathit{\theta }}_{\mathrm{f}x}}$ and $k_{{\mathit{\theta }}_{\mathrm{f}x}}$ are the damping and torsional pendulum vibration stiffness of face gear in the x axis, respectively. $L_A=\sqrt{r_{\mathrm{ca}}^{\mathrm{2}}-r_{\mathrm{cb}}^{\mathrm{2}}}-\mathit{\epsilon }{s}_{\mathrm{cb}}+n_{\mathrm{c}}{r}_{\mathrm{cb}}t$ is the friction arm. r_cb and r_ca are the base radius and addendum radius of the cylindrical gear, respectively. ε is the overlap ratio of the gear pair. s_cb is the base tooth pitch. n_c is the speed of the cylindrical gear. ${\overline{f}}_{\mathit{\mu }}=\text{sign}(L_A-r_{\mathrm{cb}}\tan {\mathit{\alpha }}_{\mathrm{n}})$ is the direction function of the frictional force.

3 Solving of coupled vibration differential equation of face gear

In order to reduce the vibrational degree of freedom in the system, it is necessary to simultaneous the torsional vibration equations:

$$\begin{array}{}\text{(6)}& \left\{\begin{array}{l}{\mathbf{J}}_{\mathrm{c}y}{\ddot{\mathit{\theta }}}_{\mathrm{c}y}={\mathbf{T}}_{\mathrm{c}}-{\mathbf{F}}_{\mathrm{n}}{r}_{\mathrm{cm}}+\mathit{\mu }{\overline{f}}_{\mathit{\mu }}{\mathbf{L}}_A{\mathbf{F}}_{\mathrm{n}}\sin {\mathit{\alpha }}_{\mathrm{n}}\\ {\mathbf{J}}_{\mathrm{f}z}{\ddot{\mathit{\theta }}}_{\mathrm{f}z}={\mathbf{F}}_{\mathrm{n}}{r}_{\mathrm{fm}}-{\mathbf{T}}_{\mathrm{f}}+\mathit{\mu }{\overline{f}}_{\mathit{\mu }}{r}_{\mathrm{fm}}{\mathbf{F}}_{\mathrm{n}}\cos {\mathit{\alpha }}_{\mathrm{n}}\end{array}\right.\end{array}$$

After rewriting:

$$\begin{array}{}\text{(7)}& \left\{\begin{array}{l}{r}_{\mathrm{cm}}{\ddot{\mathit{\theta }}}_{\mathrm{c}y}=\frac{{\mathbf{T}}_{\mathrm{c}}-{\mathbf{F}}_{\mathrm{n}}{r}_{\mathrm{cm}}+\mathit{\mu }{\overline{f}}_{\mathit{\mu }}{\mathbf{L}}_A{\mathbf{F}}_{\mathrm{n}}\sin {\mathit{\alpha }}_{\mathrm{n}}}{{\mathbf{J}}_{\mathrm{c}y}}{r}_{\mathrm{cm}}\\ r_{\mathrm{fm}}{\ddot{\mathit{\theta }}}_{\mathrm{f}z}=\frac{{\mathbf{F}}_{\mathrm{n}}{r}_{\mathrm{fm}}-{\mathbf{T}}_{\mathrm{f}}+\mathit{\mu }{\overline{f}}_{\mathit{\mu }}{r}_{\mathrm{fm}}{\mathbf{F}}_{\mathrm{n}}\cos {\mathit{\alpha }}_{\mathrm{n}}}{{\mathbf{J}}_{\mathrm{f}z}}{r}_{\mathrm{fm}}\end{array}\right.\end{array}$$

Hence,

$$\begin{array}{}\text{(8)}& \begin{array}{rl}{r}_{\mathrm{cm}}{\ddot{\mathit{\theta }}}_{\mathrm{c}y}-r_{\mathrm{fm}}{\ddot{\mathit{\theta }}}_{\mathrm{f}z}& ={\displaystyle \frac{{\mathbf{T}}_{\mathrm{c}}-{\mathbf{F}}_{\mathrm{n}}{r}_{\mathrm{cm}}}{{\mathbf{J}}_{\mathrm{c}y}}}{r}_{\mathrm{cm}}-{\displaystyle \frac{{\mathbf{F}}_{\mathrm{n}}{r}_{\mathrm{fm}}-{\mathbf{T}}_{\mathrm{f}}}{{\mathbf{J}}_{\mathrm{f}z}}}{r}_{\mathrm{fm}}\\ & +{\displaystyle \frac{\mathit{\mu }{\overline{f}}_{\mathit{\mu }}{\mathbf{L}}_A{\mathbf{F}}_{\mathrm{n}}\sin {\mathit{\alpha }}_{\mathrm{n}}}{{\mathbf{J}}_{\mathrm{c}y}}}{r}_{\mathrm{cm}}\\ & -{\displaystyle \frac{\mathit{\mu }{\overline{f}}_{\mathit{\mu }}{r}_{\mathrm{fm}}{\mathbf{F}}_{\mathrm{n}}\cos {\mathit{\alpha }}_{\mathrm{n}}}{{\mathbf{J}}_{\mathrm{f}z}}}{r}_{\mathrm{fm}}\end{array}\end{array}$$

The first and second derivative ${\dot{d}}_{\mathrm{n}}$, ${\ddot{d}}_{\mathrm{n}}$ of displacement d_n in the normal direction of mesh points are introduced:

$$\begin{array}{}\text{(9)}& \left\{\begin{array}{ll}{\dot{d}}_{\mathrm{n}}& =({\dot{x}}_{\mathrm{c}}-{\dot{x}}_{\mathrm{f}})\cos {\mathit{\alpha }}_{\mathrm{n}}+(r_{\mathrm{cm}}{\dot{\mathit{\theta }}}_{\mathrm{c}y}-r_{\mathrm{fm}}{\dot{\mathit{\theta }}}_{\mathrm{f}z})\cos {\mathit{\alpha }}_{\mathrm{n}}\\ & +({\dot{z}}_{\mathrm{c}}-{\dot{z}}_{\mathrm{f}}-r_{\mathrm{fm}}{\dot{\mathit{\theta }}}_{\mathrm{f}x})\sin {\mathit{\alpha }}_{\mathrm{n}}-{\dot{e}}_{\mathrm{n}}\left(t\right)\\ {\ddot{d}}_{\mathrm{n}}& =({\ddot{x}}_{\mathrm{c}}-{\ddot{x}}_{\mathrm{f}})\cos {\mathit{\alpha }}_{\mathrm{n}}+(r_{\mathrm{cm}}{\ddot{\mathit{\theta }}}_{\mathrm{c}y}-r_{\mathrm{fm}}{\ddot{\mathit{\theta }}}_{\mathrm{f}z})\cos {\mathit{\alpha }}_{\mathrm{n}}\\ & +({\ddot{z}}_{\mathrm{c}}-{\ddot{z}}_{\mathrm{f}}-r_{\mathrm{fm}}{\ddot{\mathit{\theta }}}_{\mathrm{f}x})\sin {\mathit{\alpha }}_{\mathrm{n}}-{\ddot{e}}_{\mathrm{n}}\left(t\right)\end{array}\right.\end{array}$$

Equation (8) is substituted and rearranged as:

$$\begin{array}{}\text{(10)}& \begin{array}{rl}{\ddot{d}}_{\mathrm{n}}& =({\ddot{x}}_{\mathrm{c}}-{\ddot{x}}_{\mathrm{f}})\cos {\mathit{\alpha }}_{\mathrm{n}}+\left({\displaystyle \frac{{\mathbf{T}}_{\mathrm{c}}-{\mathbf{F}}_{\mathrm{n}}{r}_{\mathrm{cm}}}{{\mathbf{J}}_{\mathrm{c}y}}}{r}_{\mathrm{cm}}-{\displaystyle \frac{{\mathbf{F}}_{\mathrm{n}}{r}_{\mathrm{fm}}-{\mathbf{T}}_{\mathrm{f}}}{{\mathbf{J}}_{\mathrm{f}z}}}{r}_{\mathrm{fm}}\right.\\ & \left.+{\displaystyle \frac{\mathit{\mu }{\overline{f}}_{\mathit{\mu }}{\mathbf{L}}_A{\mathbf{F}}_{\mathrm{n}}\sin {\mathit{\alpha }}_{\mathrm{n}}}{{\mathbf{J}}_{\mathrm{c}y}}}{r}_{\mathrm{cm}}-{\displaystyle \frac{\mathit{\mu }{\overline{f}}_{\mathit{\mu }}{r}_{\mathrm{fm}}{\mathbf{F}}_{\mathrm{n}}\cos {\mathit{\alpha }}_{\mathrm{n}}}{{\mathbf{J}}_{\mathrm{f}z}}}{r}_{\mathrm{fm}}\right)\cos {\mathit{\alpha }}_{\mathrm{n}}\\ & +({\ddot{z}}_{\mathrm{c}}-{\ddot{z}}_{\mathrm{f}}-r_{\mathrm{fm}}{\ddot{\mathit{\theta }}}_{\mathrm{f}x})\sin {\mathit{\alpha }}_{\mathrm{n}}-{\ddot{e}}_{\mathrm{n}}\left(t\right)\end{array}\end{array}$$

M_e is substituted and rearranged as:

$$\begin{array}{}\text{(11)}& \begin{array}{rl}{\mathbf{M}}_e{\ddot{d}}_{\mathrm{n}}& -{\mathbf{M}}_e({\ddot{x}}_{\mathrm{c}}-{\ddot{x}}_{\mathrm{f}})\cos {\mathit{\alpha }}_{\mathrm{n}}-{\mathbf{M}}_e({\ddot{z}}_{\mathrm{c}}-{\ddot{z}}_{\mathrm{f}})\sin {\mathit{\alpha }}_{\mathrm{n}}\\ & +{\mathbf{M}}_e{\ddot{e}}_{\mathrm{n}}\left(t\right)={\displaystyle \frac{{\mathbf{M}}_e{\mathbf{T}}_{\mathrm{c}}{r}_{\mathrm{cm}}\cos {\mathit{\alpha }}_{\mathrm{n}}}{{\mathbf{J}}_{\mathrm{c}y}}}+{\displaystyle \frac{{\mathbf{M}}_e{\mathbf{T}}_{\mathrm{f}}{r}_{\mathrm{fm}}\cos {\mathit{\alpha }}_{\mathrm{n}}}{{\mathbf{J}}_{\mathrm{f}z}}}\\ & +{\displaystyle \frac{{\mathbf{M}}_e\mathit{\mu }{\overline{f}}_{\mathit{\mu }}{\mathbf{L}}_A{r}_{\mathrm{cm}}{\mathbf{F}}_{\mathrm{n}}\sin {\mathit{\alpha }}_{\mathrm{n}}\cos {\mathit{\alpha }}_{\mathrm{n}}}{{\mathbf{J}}_{\mathrm{c}y}}}\\ & -{\displaystyle \frac{{\mathbf{M}}_e\mathit{\mu }{\overline{f}}_{\mathit{\mu }}{r}_{\mathrm{fm}}^{\mathrm{2}}{\mathbf{F}}_{\mathrm{n}}{\cos }^{\mathrm{2}}{\mathit{\alpha }}_{\mathrm{n}}}{{\mathbf{J}}_{\mathrm{f}z}}}\\ & -\cos {\mathit{\alpha }}_{\mathrm{n}}{\mathbf{F}}_{\mathrm{n}}-{\mathbf{M}}_e{r}_{\mathrm{fm}}{\ddot{\mathit{\theta }}}_{\mathrm{f}x}\sin {\mathit{\alpha }}_{\mathrm{n}}\end{array}\end{array}$$

Let ${\mathbf{F}}_{\mathrm{c}y}^{\prime }=\frac{{\mathbf{M}}_e{\mathbf{T}}_{\mathrm{c}}{r}_{\mathrm{cm}}\cos {\mathit{\alpha }}_{\mathrm{n}}}{{\mathbf{J}}_{\mathrm{c}y}}$, ${\mathbf{F}}_{\mathrm{f}z}^{\prime }=\frac{{\mathbf{M}}_e{\mathbf{T}}_{\mathrm{f}}{r}_{\mathrm{fm}}\cos {\mathit{\alpha }}_{\mathrm{n}}}{{\mathbf{J}}_{\mathrm{f}z}}$, ${\mathbf{F}}_{\mathrm{c}y}^{\prime \prime }=\frac{{\mathbf{M}}_e\mathit{\mu }{\overline{f}}_{\mathit{\mu }}{\mathbf{L}}_A{r}_{\mathrm{cm}}{\mathbf{F}}_{\mathrm{n}}\sin {\mathit{\alpha }}_{\mathrm{n}}\cos {\mathit{\alpha }}_{\mathrm{n}}}{{\mathbf{J}}_{\mathrm{c}y}}$, ${\mathbf{F}}_{\mathrm{f}z}^{\prime \prime }=\frac{{\mathbf{M}}_e\mathit{\mu }{\overline{f}}_{\mathit{\mu }}{r}_{\mathrm{fm}}^{\mathrm{2}}{\mathbf{F}}_{\mathrm{n}}{\cos }^{\mathrm{2}}{\mathit{\alpha }}_{\mathrm{n}}}{{\mathbf{J}}_{\mathrm{f}z}}$ be substituted and simplified as:

$$\begin{array}{}\text{(12)}& \begin{array}{rl}{\mathbf{M}}_e{\ddot{d}}_{\mathrm{n}}& -{\mathbf{M}}_e({\ddot{x}}_{\mathrm{c}}-{\ddot{x}}_{\mathrm{f}})\cos {\mathit{\alpha }}_{\mathrm{n}}-{\mathbf{M}}_e({\ddot{z}}_{\mathrm{c}}-{\ddot{z}}_{\mathrm{f}})\sin {\mathit{\alpha }}_{\mathrm{n}}+{\mathbf{M}}_e{\ddot{e}}_{\mathrm{n}}\left(t\right)\\ & ={\mathbf{F}}_{\mathrm{c}y}^{\prime }+{\mathbf{F}}_{\mathrm{f}z}^{\prime }+{\mathbf{F}}_{\mathrm{c}y}^{\prime \prime }-{\mathbf{F}}_{\mathrm{f}z}^{\prime \prime }-\cos {\mathit{\alpha }}_{\mathrm{n}}{\mathbf{F}}_{\mathrm{n}}-{\mathbf{M}}_e{r}_{\mathrm{fm}}{\ddot{\mathit{\theta }}}_{\mathrm{f}x}\sin {\mathit{\alpha }}_{\mathrm{n}}\end{array}\end{array}$$

Then, the coupled vibration equation can be reduced to six degrees of freedom:

$$\begin{array}{}\text{(13)}& \left\{\begin{array}{l}{\mathbf{M}}_{\mathrm{c}}{\ddot{x}}_{\mathrm{c}}+c_{\mathrm{c}x}{\dot{x}}_{\mathrm{c}}+k_{\mathrm{c}x}{x}_{\mathrm{c}}=-{\mathbf{F}}_x-\mathit{\mu }{\overline{f}}_{\mathit{\mu }}{\mathbf{F}}_{\mathrm{n}}\sin {\mathit{\alpha }}_{\mathrm{n}}\\ {\mathbf{M}}_{\mathrm{c}}{\ddot{z}}_{\mathrm{c}}+c_{\mathrm{c}z}{\dot{z}}_{\mathrm{c}}+k_{\mathrm{c}z}{z}_{\mathrm{c}}=-{\mathbf{F}}_z-\mathit{\mu }{\overline{f}}_{\mathit{\mu }}{\mathbf{F}}_{\mathrm{n}}\cos {\mathit{\alpha }}_{\mathrm{n}}\\ {\mathbf{M}}_{\mathrm{f}}{\ddot{x}}_{\mathrm{f}}+c_{\mathrm{f}x}{\dot{x}}_{\mathrm{f}}+k_{\mathrm{f}x}{x}_{\mathrm{f}}={\mathbf{F}}_x+\mathit{\mu }{\overline{f}}_{\mathit{\mu }}{\mathbf{F}}_{\mathrm{n}}\sin {\mathit{\alpha }}_{\mathrm{n}}\\ {\mathbf{M}}_{\mathrm{f}}{\ddot{z}}_{\mathrm{f}}+c_{\mathrm{f}z}{\dot{z}}_{\mathrm{f}}+k_{\mathrm{f}z}{z}_{\mathrm{f}}={\mathbf{F}}_z+\mathit{\mu }{\overline{f}}_{\mathit{\mu }}{\mathbf{F}}_{\mathrm{n}}\cos {\mathit{\alpha }}_{\mathrm{n}}\\ {\mathbf{J}}_{\mathrm{f}x}{\ddot{\mathit{\theta }}}_{\mathrm{f}x}+c_{{\mathit{\theta }}_{\mathrm{f}x}}{\dot{\mathit{\theta }}}_{\mathrm{f}x}+k_{{\mathit{\theta }}_{\mathrm{f}x}}{\mathit{\theta }}_{\mathrm{f}x}={\mathbf{F}}_z{r}_{\mathrm{fm}}-\mathit{\mu }{\overline{f}}_{\mathit{\mu }}{\mathbf{F}}_{\mathrm{n}}\sin {\mathit{\alpha }}_{\mathrm{n}}\\ {\mathbf{M}}_e{\ddot{d}}_{\mathrm{n}}-{\mathbf{M}}_e({\ddot{x}}_{\mathrm{c}}-{\ddot{x}}_{\mathrm{f}})\cos {\mathit{\alpha }}_{\mathrm{n}}-{\mathbf{M}}_e({\ddot{z}}_{\mathrm{c}}-{\ddot{z}}_{\mathrm{f}})\sin {\mathit{\alpha }}_{\mathrm{n}}+{\mathbf{M}}_e{\ddot{e}}_{\mathrm{n}}\left(t\right)\\ ={\mathbf{F}}_{\mathrm{c}y}^{\prime }+{\mathbf{F}}_{\mathrm{f}z}^{\prime }+{\mathbf{F}}_{\mathrm{c}y}^{\prime \prime }-{\mathbf{F}}_{\mathrm{f}z}^{\prime \prime }-\cos {\mathit{\alpha }}_{\mathrm{n}}{\mathbf{F}}_{\mathrm{n}}-{\mathbf{M}}_e{r}_{\mathrm{fm}}{\ddot{\mathit{\theta }}}_{\mathrm{f}x}\sin {\mathit{\alpha }}_{\mathrm{n}}\end{array}\right.\end{array}$$

In the process of solving the differential equation of the gear coupling vibration, the torsional pendulum vibration has a little influence on the gear transmission system, so it can be neglected. The coupled vibration equation was reduced to five degrees of freedom. Equation (3) is introduced and rearranged as:

$$\begin{array}{}\text{(14)}& \left\{\begin{array}{l}{\mathbf{M}}_{\mathrm{c}}{\ddot{x}}_{\mathrm{c}}+c_{\mathrm{c}x}{\dot{x}}_{\mathrm{c}}+k_{\mathrm{c}x}{x}_{\mathrm{c}}=-(\cos {\mathit{\alpha }}_{\mathrm{n}}+\mathit{\mu }{\overline{f}}_{\mathit{\mu }}\sin {\mathit{\alpha }}_{\mathrm{n}}){\mathbf{F}}_{\mathrm{n}}\\ {\mathbf{M}}_{\mathrm{c}}{\ddot{z}}_{\mathrm{c}}+c_{\mathrm{c}z}{\dot{z}}_{\mathrm{c}}+k_{\mathrm{c}z}{z}_{\mathrm{c}}=-(\sin {\mathit{\alpha }}_{\mathrm{n}}-\mathit{\mu }{\overline{f}}_{\mathit{\mu }}\cos {\mathit{\alpha }}_{\mathrm{n}}){\mathbf{F}}_{\mathrm{n}}\\ {\mathbf{M}}_{\mathrm{f}}{\ddot{x}}_{\mathrm{f}}+c_{\mathrm{f}x}{\dot{x}}_{\mathrm{f}}+k_{\mathrm{f}x}{x}_{\mathrm{f}}=(\cos {\mathit{\alpha }}_{\mathrm{n}}+\mathit{\mu }{\overline{f}}_{\mathit{\mu }}\sin {\mathit{\alpha }}_{\mathrm{n}}){\mathbf{F}}_{\mathrm{n}}\\ {\mathbf{M}}_{\mathrm{f}}{\ddot{z}}_{\mathrm{f}}+c_{\mathrm{f}z}{\dot{z}}_{\mathrm{f}}+k_{\mathrm{f}z}{z}_{\mathrm{f}}=(\sin {\mathit{\alpha }}_{\mathrm{n}}+\mathit{\mu }{\overline{f}}_{\mathit{\mu }}\cos {\mathit{\alpha }}_{\mathrm{n}}){\mathbf{F}}_{\mathrm{n}}\\ {\mathbf{M}}_e{\ddot{d}}_{\mathrm{n}}-{\mathbf{M}}_e({\ddot{x}}_{\mathrm{c}}-{\ddot{x}}_{\mathrm{f}})\cos {\mathit{\alpha }}_{\mathrm{n}}-{\mathbf{M}}_e({\ddot{z}}_{\mathrm{c}}-{\ddot{z}}_{\mathrm{f}})\sin {\mathit{\alpha }}_{\mathrm{n}}\\ +{\mathbf{M}}_e{\ddot{e}}_{\mathrm{n}}\left(t\right)={\mathbf{F}}_{\mathrm{c}y}^{\prime }+{\mathbf{F}}_{\mathrm{f}z}^{\prime }+{\mathbf{F}}_{\mathrm{c}y}^{\prime \prime }-{\mathbf{F}}_{\mathrm{f}z}^{\prime \prime }-\cos {\mathit{\alpha }}_{\mathrm{n}}{\mathbf{F}}_{\mathrm{n}}\end{array}\right.\end{array}$$

Due to the large magnitude of differences between the physical quantities, a numerical solution is extremely difficult, and it is also challenging to select the step size and control the error. Therefore, it is necessary to adopt the dimensionless normalization to the coupled vibration equation. Assuming that the natural frequency of the dynamic model of face gear pair is:

$$\begin{array}{}\text{(15)}& N_{\mathrm{f}}=\sqrt{{\overline{k}}_{\mathrm{cm}}/{\mathbf{M}}_e}\end{array}$$

The dimensionless excitation frequency N_D and the time independent variable τ are respectively defined as:

$$\begin{array}{}\text{(16)}& \left\{\begin{array}{l}{N}_{\mathrm{D}}={\mathit{\omega }}_i/N_{\mathrm{f}}\\ \mathit{\tau }=N_{\mathrm{f}}t\end{array}\right.\end{array}$$

The results of the dimensionless transformation of the new 5 degrees of freedom coupled vibration equation are:

$$\begin{array}{}\text{(17)}& \left\{\begin{array}{l}\begin{array}{ll}\frac{d^{\mathrm{2}}{x}_{\mathrm{dc}}}{\mathrm{d}{\mathit{\tau }}^{\mathrm{2}}}& +\frac{c_{\mathrm{c}x}}{{\mathbf{M}}_{\mathrm{dc}}{N}_{\mathrm{f}}}\frac{\mathrm{d}{x}_{\mathrm{dc}}}{\mathrm{d}\mathit{\tau }}+\frac{k_{\mathrm{c}x}}{{\mathbf{M}}_{\mathrm{dc}}{N}_{\mathrm{f}}^{\mathrm{2}}}{x}_{\mathrm{dc}}\\ & -(\cos {\mathit{\alpha }}_{\mathrm{n}}+\mathit{\mu }{\overline{f}}_{\mathit{\mu }}\sin {\mathit{\alpha }}_{\mathrm{n}})\frac{k_{\mathrm{n}f}\left(\mathit{\tau }\right)}{{\mathbf{M}}_{\mathrm{dc}}{N}_{\mathrm{f}}^{\mathrm{2}}}{f}_{\mathrm{b}}\left(d_{\mathrm{d}}\right)\\ & -(\cos {\mathit{\alpha }}_{\mathrm{n}}+\mathit{\mu }{\overline{f}}_{\mathit{\mu }}\sin {\mathit{\alpha }}_{\mathrm{n}})\frac{c_{\mathrm{cf}}}{{\mathbf{M}}_{\mathrm{dc}}{N}_{\mathrm{f}}}\frac{d{d}_{\mathrm{d}}}{\mathrm{d}\mathit{\tau }}=\mathrm{0}\end{array}\\ \begin{array}{ll}\frac{d^{\mathrm{2}}{z}_{\mathrm{dc}}}{\mathrm{d}{\mathit{\tau }}^{\mathrm{2}}}& +\frac{c_{\mathrm{c}z}}{{\mathbf{M}}_{\mathrm{dc}}{N}_{\mathrm{f}}}\frac{\mathrm{d}{z}_{\mathrm{dc}}}{\mathrm{d}\mathit{\tau }}+\frac{k_{\mathrm{c}z}}{{\mathbf{M}}_{\mathrm{dc}}{N}_{\mathrm{f}}^{\mathrm{2}}}{z}_{\mathrm{dc}}\\ & +(\sin {\mathit{\alpha }}_{\mathrm{n}}-\mathit{\mu }{\overline{f}}_{\mathit{\mu }}\cos {\mathit{\alpha }}_{\mathrm{n}})\frac{k_{\mathrm{n}f}\left(\mathit{\tau }\right)}{{\mathbf{M}}_{\mathrm{dc}}{N}_{\mathrm{f}}^{\mathrm{2}}}{f}_{\mathrm{b}}\left(d_{\mathrm{d}}\right)\\ & +(\sin {\mathit{\alpha }}_{\mathrm{n}}-\mathit{\mu }{\overline{f}}_{\mathit{\mu }}\cos {\mathit{\alpha }}_{\mathrm{n}})\frac{c_{\mathrm{c}f}}{{\mathbf{M}}_{\mathrm{dc}}{N}_{\mathrm{f}}}\frac{\mathrm{d}{d}_{\mathrm{d}}}{\mathrm{d}\mathit{\tau }}=\mathrm{0}\end{array}\\ \begin{array}{ll}\frac{d^{\mathrm{2}}{x}_{\mathrm{df}}}{\mathrm{d}{\mathit{\tau }}^{\mathrm{2}}}& +\frac{c_{\mathrm{f}x}}{{\mathbf{M}}_{\mathrm{df}}{N}_{\mathrm{f}}}\frac{\mathrm{d}{x}_{\mathrm{df}}}{\mathrm{d}\mathit{\tau }}+\frac{k_{\mathrm{f}x}}{{\mathbf{M}}_{\mathrm{df}}{N}_{\mathrm{f}}^{\mathrm{2}}}{x}_{\mathrm{df}}\\ & -(\cos {\mathit{\alpha }}_{\mathrm{n}}+\mathit{\mu }{\overline{f}}_{\mathit{\mu }}\sin {\mathit{\alpha }}_{\mathrm{n}})\frac{k_{\mathrm{n}f}\left(\mathit{\tau }\right)}{{\mathbf{M}}_{\mathrm{df}}{N}_{\mathrm{f}}^{\mathrm{2}}}{f}_{\mathrm{b}}\left(d_{\mathrm{d}}\right)\\ & -(\cos {\mathit{\alpha }}_{\mathrm{n}}+\mathit{\mu }{\overline{f}}_{\mathit{\mu }}\sin {\mathit{\alpha }}_{\mathrm{n}})\frac{c_{\mathrm{cf}}}{{\mathbf{M}}_{\mathrm{df}}{N}_{\mathrm{f}}}\frac{\mathrm{d}{d}_{\mathrm{d}}}{\mathrm{d}\mathit{\tau }}=\mathrm{0}\end{array}\\ \begin{array}{ll}\frac{d^{\mathrm{2}}{z}_{\mathrm{df}}}{\mathrm{d}{\mathit{\tau }}^{\mathrm{2}}}& +\frac{c_{\mathrm{f}z}}{{\mathbf{M}}_{\mathrm{df}}{N}_{\mathrm{f}}}\frac{\mathrm{d}{z}_{\mathrm{df}}}{\mathrm{d}\mathit{\tau }}+\frac{k_{\mathrm{f}z}}{{\mathbf{M}}_{\mathrm{df}}{N}_{\mathrm{f}}^{\mathrm{2}}}{z}_{\mathrm{df}}\\ & +(\sin {\mathit{\alpha }}_{\mathrm{n}}-\mathit{\mu }{\overline{f}}_{\mathit{\mu }}\cos {\mathit{\alpha }}_{\mathrm{n}})\frac{k_{\mathrm{nf}}\left(\mathit{\tau }\right)}{{\mathbf{M}}_{\mathrm{df}}{N}_{\mathrm{f}}^{\mathrm{2}}}{f}_{\mathrm{b}}\left(d_{\mathrm{d}}\right)\\ & +(\sin {\mathit{\alpha }}_{\mathrm{n}}-\mathit{\mu }{\overline{f}}_{\mathit{\mu }}\cos {\mathit{\alpha }}_{\mathrm{n}})\frac{c_{\mathrm{cf}}}{{\mathbf{M}}_{\mathrm{df}}{N}_{\mathrm{f}}}\frac{\mathrm{d}{d}_{\mathrm{d}}}{\mathrm{d}\mathit{\tau }}=\mathrm{0}\end{array}\\ \begin{array}{ll}\frac{d^{\mathrm{2}}{d}_{\mathrm{d}}}{\mathrm{d}{\mathit{\tau }}^{\mathrm{2}}}& -\left(\frac{d^{\mathrm{2}}{x}_{\mathrm{dc}}}{\mathrm{d}{\mathit{\tau }}^{\mathrm{2}}}-\frac{d^{\mathrm{2}}{x}_{\mathrm{df}}}{\mathrm{d}{\mathit{\tau }}^{\mathrm{2}}}\right)\cos {\mathit{\alpha }}_{\mathrm{n}}-\left(\frac{d^{\mathrm{2}}{z}_{\mathrm{dc}}}{\mathrm{d}{\mathit{\tau }}^{\mathrm{2}}}-\frac{d^{\mathrm{2}}{z}_{\mathrm{df}}}{\mathrm{d}{\mathit{\tau }}^{\mathrm{2}}}\right)\sin {\mathit{\alpha }}_{\mathrm{n}}\\ & =\frac{{\mathbf{F}}_{\mathrm{c}y}^{\prime }+{\mathbf{F}}_{\mathrm{f}z}^{\prime }+{\mathbf{F}}_{\mathrm{c}y}^{\prime \prime }-{\mathbf{F}}_{\mathrm{f}z}^{\prime \prime }}{{\mathbf{M}}_{\mathrm{d}e}{N}_{\mathrm{f}}^{\mathrm{2}}}-\frac{k_{\mathrm{nf}}\left(\mathit{\tau }\right)\cos {\mathit{\alpha }}_{\mathrm{n}}}{{\mathbf{M}}_{\mathrm{d}e}{N}_{\mathrm{f}}^{\mathrm{2}}}{f}_{\mathrm{b}}\left(d_{\mathrm{d}}\right)\\ & -\frac{c_{\mathrm{cf}}\left(\mathit{\tau }\right)\cos {\mathit{\alpha }}_{\mathrm{n}}}{{\mathbf{M}}_{\mathrm{d}e}{N}_{\mathrm{f}}}\frac{\mathrm{d}{d}_{\mathrm{d}}}{\mathrm{d}\mathit{\tau }}-\frac{\mathrm{1}}{N_{\mathrm{f}}^{\mathrm{2}}}{\ddot{e}}_{\mathrm{n}f}\left(\mathit{\tau }\right)\end{array}\end{array}\right.\end{array}$$

The static integrated transmission error equation is introduced and derivative is taken:

$$\begin{array}{}\text{(18)}& \left\{\begin{array}{l}{\dot{e}}_{\mathrm{nf}}\left(\mathit{\tau }\right)=e_r{\mathit{\omega }}_i\cos ({\mathit{\omega }}_i\mathit{\tau }+{\mathit{\phi }}_{ipa})\\ {\ddot{e}}_{\mathrm{n}f}\left(\mathit{\tau }\right)=-e_r{\mathit{\omega }}_i^{\mathrm{2}}\sin ({\mathit{\omega }}_i\mathit{\tau }+{\mathit{\phi }}_{ipa})\end{array}\right.\end{array}$$

Let $A_{\mathrm{c}x}=\frac{c_{\mathrm{c}x}}{\mathrm{2}{\mathbf{M}}_{\mathrm{dc}}{N}_{\mathrm{f}}}$, $A_{\mathrm{dc}}=\frac{c_{\mathrm{cf}}}{\mathrm{2}{\mathbf{M}}_{\mathrm{dc}}{N}_{\mathrm{f}}}$, $B_{\mathrm{c}x}=\frac{k_{\mathrm{c}x}}{{\mathbf{M}}_{\mathrm{dc}}{N}_{\mathrm{f}}^{\mathrm{2}}}$, $B_{\mathrm{dc}}=\frac{k_{\mathrm{nf}}\left(\mathit{\tau }\right)}{{\mathbf{M}}_{\mathrm{dc}}{N}_{\mathrm{f}}^{\mathrm{2}}}$; $A_{\mathrm{c}z}=\frac{c_{\mathrm{c}z}}{\mathrm{2}{\mathbf{M}}_{\mathrm{dc}}{N}_{\mathrm{f}}}$, $B_{\mathrm{c}z}=\frac{k_{\mathrm{c}z}}{{\mathbf{M}}_{\mathrm{dc}}{N}_{\mathrm{f}}^{\mathrm{2}}}$; $A_{\mathrm{f}x}=\frac{c_{\mathrm{f}x}}{\mathrm{2}{\mathbf{M}}_{\mathrm{df}}{N}_{\mathrm{f}}}$, $A_{\mathrm{df}}=\frac{c_{\mathrm{cf}}}{\mathrm{2}{\mathbf{M}}_{\mathrm{df}}{N}_{\mathrm{f}}}$, $B_{\mathrm{f}x}=\frac{k_{\mathrm{f}x}}{{\mathbf{M}}_{\mathrm{df}}{N}_{\mathrm{f}}^{\mathrm{2}}}$, $B_{\mathrm{df}}=\frac{k_{\mathrm{n}f}\left(\mathit{\tau }\right)}{{\mathbf{M}}_{\mathrm{df}}{N}_{\mathrm{f}}^{\mathrm{2}}}$; $A_{\mathrm{f}z}=\frac{c_{\mathrm{f}z}}{\mathrm{2}{\mathbf{M}}_{\mathrm{df}}{N}_{\mathrm{f}}}$, $B_{\mathrm{f}z}=\frac{k_{\mathrm{f}z}}{{\mathbf{M}}_{\mathrm{df}}{N}_{\mathrm{f}}^{\mathrm{2}}}$; $B_{\mathrm{d}e}=\frac{k_{\mathrm{nf}}\left(\mathit{\tau }\right)}{{\mathbf{M}}_{\mathrm{d}e}{N}_{\mathrm{f}}^{\mathrm{2}}}$, $A_{\mathrm{d}e}=\frac{c_{\mathrm{cf}}\left(\mathit{\tau }\right)}{\mathrm{2}{\mathbf{M}}_{\mathrm{d}e}{N}_{\mathrm{f}}}$, $C=\frac{{\mathbf{F}}_{\mathrm{c}y}^{\prime }+{\mathbf{F}}_{\mathrm{f}z}^{\prime }+{\mathbf{F}}_{\mathrm{c}y}^{\prime \prime }-{\mathbf{F}}_{\mathrm{f}z}^{\prime \prime }}{{\mathbf{M}}_{\mathrm{d}e}{N}_{\mathrm{f}}^{\mathrm{2}}}$, $D=\frac{e_r{\mathit{\omega }}_i^{\mathrm{2}}\sin ({\mathit{\omega }}_i\mathit{\tau }+{\mathit{\phi }}_{ipa})}{b_k}$.

Then, the dimensionless coupling vibration differential equation of the face gear can be expressed as:

$$\begin{array}{}\text{(19)}& \left\{\begin{array}{ll}\frac{d^{\mathrm{2}}{x}_{\mathrm{dc}}}{\mathrm{d}{\mathit{\tau }}^{\mathrm{2}}}& +\mathrm{2}{A}_{\mathrm{c}x}\frac{\mathrm{d}{x}_{\mathrm{dc}}}{\mathrm{d}\mathit{\tau }}+B_{\mathrm{c}x}{x}_{\mathrm{dc}}\\ & -(\cos {\mathit{\alpha }}_{\mathrm{n}}+\mathit{\mu }{\overline{f}}_{\mathit{\mu }}\sin {\mathit{\alpha }}_{\mathrm{n}})B_{\mathrm{dc}}{f}_{\mathrm{b}}\left(d_{\mathrm{d}}\right)\\ & -\mathrm{2}(\cos {\mathit{\alpha }}_{\mathrm{n}}+\mathit{\mu }{\overline{f}}_{\mathit{\mu }}\sin {\mathit{\alpha }}_{\mathrm{n}})A_{\mathrm{dc}}\frac{\mathrm{d}{d}_{\mathrm{d}}}{\mathrm{d}\mathit{\tau }}=\mathrm{0}\\ \frac{d^{\mathrm{2}}{z}_{\mathrm{dc}}}{\mathrm{d}{\mathit{\tau }}^{\mathrm{2}}}& +\mathrm{2}{A}_{\mathrm{c}z}\frac{\mathrm{d}{z}_{\mathrm{dc}}}{\mathrm{d}\mathit{\tau }}+B_{\mathrm{c}z}{z}_{\mathrm{dc}}\\ & +(\sin {\mathit{\alpha }}_{\mathrm{n}}-\mathit{\mu }{\overline{f}}_{\mathit{\mu }}\cos {\mathit{\alpha }}_{\mathrm{n}})B_{\mathrm{dc}}{f}_{\mathrm{b}}\left(d_{\mathrm{d}}\right)\\ & +\mathrm{2}(\sin {\mathit{\alpha }}_{\mathrm{n}}-\mathit{\mu }{\overline{f}}_{\mathit{\mu }}\cos {\mathit{\alpha }}_{\mathrm{n}})A_{\mathrm{dc}}\frac{\mathrm{d}{d}_{\mathrm{d}}}{\mathrm{d}\mathit{\tau }}=\mathrm{0}\\ \frac{d^{\mathrm{2}}{x}_{\mathrm{df}}}{\mathrm{d}{\mathit{\tau }}^{\mathrm{2}}}& +\mathrm{2}{A}_{\mathrm{f}x}\frac{\mathrm{d}{x}_{\mathrm{df}}}{\mathrm{d}\mathit{\tau }}+B_{\mathrm{f}x}{x}_{\mathrm{df}}\\ & -(\cos {\mathit{\alpha }}_{\mathrm{n}}+\mathit{\mu }{\overline{f}}_{\mathit{\mu }}\sin {\mathit{\alpha }}_{\mathrm{n}})B_{\mathrm{df}}{f}_{\mathrm{b}}\left(d_{\mathrm{d}}\right)\\ & -\mathrm{2}(\cos {\mathit{\alpha }}_{\mathrm{n}}+\mathit{\mu }{\overline{f}}_{\mathit{\mu }}\sin {\mathit{\alpha }}_{\mathrm{n}})A_{\mathrm{df}}\frac{\mathrm{d}{d}_{\mathrm{d}}}{\mathrm{d}\mathit{\tau }}=\mathrm{0}\\ \frac{d^{\mathrm{2}}{z}_{\mathrm{df}}}{\mathrm{d}{\mathit{\tau }}^{\mathrm{2}}}& +\mathrm{2}{A}_{\mathrm{f}z}\frac{\mathrm{d}{z}_{\mathrm{df}}}{\mathrm{d}\mathit{\tau }}+B_{\mathrm{f}z}{z}_{\mathrm{df}}\\ & +(\sin {\mathit{\alpha }}_{\mathrm{n}}-\mathit{\mu }{\overline{f}}_{\mathit{\mu }}\cos {\mathit{\alpha }}_{\mathrm{n}})B_{\mathrm{df}}{f}_{\mathrm{b}}\left(d_{\mathrm{d}}\right)\\ & +\mathrm{2}(\sin {\mathit{\alpha }}_{\mathrm{n}}-\mathit{\mu }{\overline{f}}_{\mathit{\mu }}\cos {\mathit{\alpha }}_{\mathrm{n}})A_{\mathrm{df}}\frac{\mathrm{d}{d}_{\mathrm{d}}}{\mathrm{d}\mathit{\tau }}=\mathrm{0}\\ \frac{d^{\mathrm{2}}{d}_{\mathrm{d}}}{\mathrm{d}{\mathit{\tau }}^{\mathrm{2}}}& -\frac{d^{\mathrm{2}}{x}_{\mathrm{dc}}}{\mathrm{d}{\mathit{\tau }}^{\mathrm{2}}}\cos {\mathit{\alpha }}_{\mathrm{n}}+\frac{d^{\mathrm{2}}{x}_{\mathrm{df}}}{\mathrm{d}{\mathit{\tau }}^{\mathrm{2}}}\cos {\mathit{\alpha }}_{\mathrm{n}}-\frac{d^{\mathrm{2}}{z}_{\mathrm{dc}}}{\mathrm{d}{\mathit{\tau }}^{\mathrm{2}}}\sin {\mathit{\alpha }}_{\mathrm{n}}\\ & +\frac{d^{\mathrm{2}}{z}_{\mathrm{df}}}{\mathrm{d}{\mathit{\tau }}^{\mathrm{2}}}\sin {\mathit{\alpha }}_{\mathrm{n}}+B_{\mathrm{d}e}\cos {\mathit{\alpha }}_{\mathrm{n}}{f}_{\mathrm{b}}\left(d_{\mathrm{d}}\right)\\ & +\mathrm{2}{A}_{\mathrm{d}e}\cos {\mathit{\alpha }}_{\mathrm{n}}\frac{\mathrm{d}{d}_{\mathrm{d}}}{\mathrm{d}\mathit{\tau }}=C+D\end{array}\right.\end{array}$$

The solution of the coupled vibration differential equation of face gear usually discretizes the time. A certain step length should be selected to obtain the relationship of acceleration, velocity, and displacement at every moment. Thereafter, integration is performed directly in steps to finally obtain both the transient and steady state solutions. When the effects of damping force, alternating excitation, and gear backlash are not considered, the integral initial values of steady state vibration can be obtained using the variable step length method. However, the face gear with double crown surface has a complex tooth surface structure, which is affected by the contact characteristics of the crown structure and the gear backlash. The meshing stiffness of the point contact is time-varying, wherein, it is difficult to carry out the complete numerical analysis. Therefore, the vibration of the face gear has strong nonlinear characteristics. Even if the linear model can well approximate the vibration response of the real system, it can only be used as a simplified model that reflects the real system. It is demanding to maintain the stability and reliability under the influence of various nonlinear errors, and its vibration characteristics are also difficult to predict.

A reduced order processing for differential equations was performed by introducing the state variable λ(τ), as a result, the Eq. (19) can be transformed into 10 first-order differential equations.

$$\begin{array}{}\text{(20)}& \begin{array}{ll}\mathit{\lambda }\left(\mathit{\tau }\right)& =\mathit{\{}{\mathit{\lambda }}_{\mathrm{1}},{\mathit{\lambda }}_{\mathrm{2}},{\mathit{\lambda }}_{\mathrm{3}},{\mathit{\lambda }}_{\mathrm{4}},{\mathit{\lambda }}_{\mathrm{5}},{\mathit{\lambda }}_{\mathrm{6}},{\mathit{\lambda }}_{\mathrm{7}},{\mathit{\lambda }}_{\mathrm{8}},{\mathit{\lambda }}_{\mathrm{9}},{\mathit{\lambda }}_{\mathrm{10}}{\mathit{\}}}^{\mathrm{T}}\\ & =\mathit{\{}{x}_{\mathrm{dc}},{\dot{x}}_{\mathrm{dc}},z_{\mathrm{dc}},{\dot{z}}_{\mathrm{dc}},x_{\mathrm{df}},{\dot{x}}_{\mathrm{df}},z_{\mathrm{df}},{\dot{z}}_{\mathrm{df}},d_{\mathrm{d}},{\dot{d}}_{\mathrm{d}}{\mathit{\}}}^{\mathrm{T}}\end{array}\end{array}$$

which are:

$$\begin{array}{}\text{(21)}& \left\{\begin{array}{ll}{\dot{\mathit{\lambda }}}_{\mathrm{1}}& ={\mathit{\lambda }}_{\mathrm{2}}\\ {\dot{\mathit{\lambda }}}_{\mathrm{2}}& =-\mathrm{2}{A}_{\mathrm{c}x}{\mathit{\lambda }}_{\mathrm{2}}-B_{\mathrm{c}x}{\mathit{\lambda }}_{\mathrm{1}}+(\cos {\mathit{\alpha }}_{\mathrm{n}}\\ & +\mathit{\mu }{\overline{f}}_{\mathit{\mu }}\sin {\mathit{\alpha }}_{\mathrm{n}}\left)B_{\mathrm{dc}}{f}_{\mathrm{b}}\right({\mathit{\lambda }}_{\mathrm{9}})\\ & +\mathrm{2}(\cos {\mathit{\alpha }}_{\mathrm{n}}+\mathit{\mu }{\overline{f}}_{\mathit{\mu }}\sin {\mathit{\alpha }}_{\mathrm{n}})A_{\mathrm{dc}}{\mathit{\lambda }}_{\mathrm{10}}\\ {\dot{\mathit{\lambda }}}_{\mathrm{3}}& ={\mathit{\lambda }}_{\mathrm{4}}\\ {\dot{\mathit{\lambda }}}_{\mathrm{4}}& =-\mathrm{2}{A}_{\mathrm{c}z}{\mathit{\lambda }}_{\mathrm{4}}-B_{\mathrm{c}z}{\mathit{\lambda }}_{\mathrm{3}}\\ & -(\sin {\mathit{\alpha }}_{\mathrm{n}}-\mathit{\mu }{\overline{f}}_{\mathit{\mu }}\cos {\mathit{\alpha }}_{\mathrm{n}})B_{\mathrm{dc}}{f}_{\mathrm{b}}\left({\mathit{\lambda }}_{\mathrm{9}}\right)\\ & -\mathrm{2}(\sin {\mathit{\alpha }}_{\mathrm{n}}-\mathit{\mu }{\overline{f}}_{\mathit{\mu }}\cos {\mathit{\alpha }}_{\mathrm{n}})A_{\mathrm{dc}}{\mathit{\lambda }}_{\mathrm{10}}\\ {\dot{\mathit{\lambda }}}_{\mathrm{5}}& ={\mathit{\lambda }}_{\mathrm{6}}\\ {\dot{\mathit{\lambda }}}_{\mathrm{6}}& =-\mathrm{2}{A}_{\mathrm{f}x}{\mathit{\lambda }}_{\mathrm{6}}-B_{\mathrm{f}x}{\mathit{\lambda }}_{\mathrm{5}}\\ & +(\cos {\mathit{\alpha }}_{\mathrm{n}}+\mathit{\mu }{\overline{f}}_{\mathit{\mu }}\sin {\mathit{\alpha }}_{\mathrm{n}})B_{\mathrm{df}}{f}_{\mathrm{b}}\left({\mathit{\lambda }}_{\mathrm{9}}\right)\\ & +\mathrm{2}(\cos {\mathit{\alpha }}_{\mathrm{n}}+\mathit{\mu }{\overline{f}}_{\mathit{\mu }}\sin {\mathit{\alpha }}_{\mathrm{n}})A_{\mathrm{df}}{\mathit{\lambda }}_{\mathrm{10}}\\ {\dot{\mathit{\lambda }}}_{\mathrm{7}}& ={\mathit{\lambda }}_{\mathrm{8}}\\ {\dot{\mathit{\lambda }}}_{\mathrm{8}}& =-\mathrm{2}{A}_{\mathrm{f}z}{\mathit{\lambda }}_{\mathrm{8}}-B_{\mathrm{f}z}{\mathit{\lambda }}_{\mathrm{7}}-(\sin {\mathit{\alpha }}_{\mathrm{n}}\\ & -\mathit{\mu }{\overline{f}}_{\mathit{\mu }}\cos {\mathit{\alpha }}_{\mathrm{n}}\left)B_{\mathrm{df}}{f}_{\mathrm{b}}\right({\mathit{\lambda }}_{\mathrm{9}})\\ & -\mathrm{2}(\sin {\mathit{\alpha }}_{\mathrm{n}}-\mathit{\mu }{\overline{f}}_{\mathit{\mu }}\cos {\mathit{\alpha }}_{\mathrm{n}})A_{\mathrm{df}}{\mathit{\lambda }}_{\mathrm{10}}\\ {\dot{\mathit{\lambda }}}_{\mathrm{9}}& ={\mathit{\lambda }}_{\mathrm{10}}\\ {\dot{\mathit{\lambda }}}_{\mathrm{10}}& =\cos {\mathit{\alpha }}_{\mathrm{n}}{\dot{\mathit{\lambda }}}_{\mathrm{2}}+\sin {\mathit{\alpha }}_{\mathrm{n}}{\dot{\mathit{\lambda }}}_{\mathrm{4}}\\ & -\cos {\mathit{\alpha }}_{\mathrm{n}}{\dot{\mathit{\lambda }}}_{\mathrm{6}}-\sin {\mathit{\alpha }}_{\mathrm{n}}{\dot{\mathit{\lambda }}}_{\mathrm{8}}-B_{\mathrm{d}e}\cos {\mathit{\alpha }}_{\mathrm{n}}{f}_{\mathrm{b}}\left({\mathit{\lambda }}_{\mathrm{9}}\right)\\ & -\mathrm{2}{A}_{\mathrm{d}e}\cos {\mathit{\alpha }}_{\mathrm{n}}{\mathit{\lambda }}_{\mathrm{10}}-(C+D)\end{array}\right.\end{array}$$

4 Analysis of linear dynamic characteristics of face gear

Neglecting the time-variation of the meshing stiffness, setting as average meshing stiffness and setting the meshing backlash to 0 mm, hence, the dimensionless coupled vibration equation can be defined as a linear differential equation. The linear characteristics of face gear were analyzed by the 4–5-order Runge–Kutta algorithm, which is the ODE45 variable step integral method in MATLAB. The main structural parameters of the gears are: tooth numbers Z_c=17, $Z_{\mathrm{f}}=\mathrm{56}/\mathrm{68}/\mathrm{80}$. Modulus m=2 mm. Pressure Angle α_n=20^∘. Input torque T_c=100 N m. The dimensionless parameters of the vibration differential equation are assumed as: $A_{\mathrm{c}x}=A_{\mathrm{c}z}=\mathrm{0.0261}$, A_dc=0.093, $B_{\mathrm{c}x}=B_{\mathrm{c}z}=\mathrm{0.0809}$, $B_{\mathrm{f}x}=B_{\mathrm{f}z}=\mathrm{0.0071}$, B_dc=0.5509, $A_{\mathrm{f}x}=A_{\mathrm{f}z}=\mathrm{3.2383}\times {\mathrm{10}}^{-\mathrm{4}}$, $A_{\mathrm{df}}=\mathrm{3.8986}\times {\mathrm{10}}^{-\mathrm{4}}$, B_df=0.0231, B_de=0.8806, A_de=0.0149, $C=\mathrm{3.3982}\times {\mathrm{10}}^{-\mathrm{5}}$, and $D=\mathrm{5.1342}\times {\mathrm{10}}^{\mathrm{6}}\times \sin \left(w_i\mathit{\tau }\right)$. At the same time, the friction coefficient is defined as 0 and the load is 0. Under these conditions, the phase diagram and Poincare section were introduced to analyze the linear dynamic characteristics, as shown in Fig. 3.

Figure 3 Displacement, phase, and Poincare diagrams of linear system.

Under the condition of without load, the time-history response of the system was a simple harmonic motion. The vibration displacement of the cylindrical gear was large at startup, and tended to be stable after a certain period of time. The face gear also presented a relatively stable period simple harmonic. The phase diagram is a plurality of closed elliptic curves that is expressed as a period response. The Poincare diagram is a number of discrete points.

4.1 Influence of variable speed stage number on the linear dynamics characteristics of the system

Other parameters were unchanged, the teeth number of the face gear were switched to the 1st, 2nd or 3rd stage, respectively, the corresponding equivalent mass, rotational inertia and gear ratio change accordingly. The displacement curves, phase diagrams, and Poincare diagrams of the linear dynamics are shown in Figs. 4, 5, and 6.

Figure 4 Displacement, phase, and Poincare diagrams of 1st stage face gear pair.

Figure 5 Displacement, phase, and Poincare diagrams of 2nd stage face gear pair.

Figure 6 Displacement, phase, and Poincare diagrams of 3rd stage face gear pair.

As per Figs. 4, 5, and 6, under the conditions without external excitation, the phase diagrams of face gears in all the stages were closed curves, which indicate that all the stages had stable vibration period characteristics. After inputting torque and speed, the system reached a stable state in an extremely short time-history, and its vibration displacement was also weak, which indicates that the multistage face gears system can achieve stable meshing transmission.

4.2 Influence of friction coefficient on linear dynamic characteristics of the system

Taking the 2nd face gear pair as the object, the other parameters remained unchanged, while the friction coefficients were set to 0.1, 0.5, and 1, respectively. Inspite of free idling of the other two stages face gears, there was still friction damping. It was necessary to define additional friction coefficients in the 1st stage and 3rd stage to be 0.01, and the 2nd stage to be 0.02. The displacement curves, phase diagrams, and Poincare diagrams of the linear dynamics are shown in Figs. 7, 8, and 9.

Figure 7 Displacement, phase, and Poincare diagrams of friction coefficient 0.1.

Figure 8 Displacement, phase, and Poincare diagrams of friction coefficient 0.5.

Figure 9 Displacement, phase, and Poincare diagrams of friction coefficient 1.

As per Figs. 7, 8, and 9, the elastic deformation was generated in gears during the meshing process, which led to the collision of the newly meshed gear teeth, thus, generating the pulse power. The meshing frequency includes not only the fundamental frequency vibration but also high harmonic vibration. The friction between two meshing teeth would also induce self-excited vibration under certain conditions. The meshing process is a comprehensive transient excitation, which causes the gear to attenuate free vibration. As the friction coefficient becomes larger, the vibration displacement increasingly deviates from the meshing trajectory. When the friction coefficient was less than 0.1, the vibration characteristics changed from single cycle to double cycle, and then, to quasi-periodic cycle. When the same was greater than 0.1, the vibration displacements in the x and z directions were gradually towards the same direction under the influence of friction damping. The Poincare diagram is presented as a complex set of points. Therefore, reasonable gear materials and appropriate tooth surface processing technology should be used to enhance the surface quality and reduce the tooth surface friction.

4.3 Influence of integrated meshing stiffness on linear dynamic characteristics of the system

Taking the 2nd face gear pair as the object, other parameters remained unchanged. The friction coefficient was set to 0.01, when the input torque is set to 100, 500, and 1000 N m⁻¹ respectively. The integrated meshing stiffness of the corresponding face gears pairs were 9.852×10⁷ N m⁻¹, 5.7617×10⁷ N m⁻¹, and 4.5731×10⁷ N m⁻¹. The displacement curves, phase diagrams, and Poincare diagrams of the linear dynamics are shown in Figs. 10, 11 and 12.

Figure 10 Displacement, phase, and Poincare diagrams of comprehensive meshing stiffness 9.852×⁷ N m⁻¹.

Figure 11 Displacement, phase, and Poincare diagrams of comprehensive meshing stiffness 5.7617×10⁷ N m⁻¹.

Figure 12 Displacement, phase, and Poincare diagrams of comprehensive meshing stiffness 4.5731×10⁷ N m⁻¹.

The integrated meshing stiffness is directly related to the input torque and the load deformation. As per Figs. 10, 11 and 12, the vibration displacement amplitude increased with the reduction of the integrated meshing stiffness. Owing to the frictional damping, the degree of vibration displacement increasingly deviated from the meshing trajectory. The phase diagram is multiple spiral curves that is expressed as a quasi-periodic response, while the Poincare diagram is a number of irregular discrete points.

4.4 Influence of loading speed on linear dynamic characteristics of the system

Taking the 2nd face gear pair as the object, the friction coefficient was 0.01, when the input torque was 100 N m⁻¹, other parameters remained unchanged, the loading speed was set to 300, 500, 1000, and 3000 r min⁻¹, respectively, and the corresponding meshing angle frequency changed accordingly.

After inputting speed excitation, the linear system generated periodic oscillations, which then decayed to a stable vibration state. Larger the input speed, longer was the oscillation period, and larger was the vibration displacement amplitude. The phase diagram presents stable multiple spiral curves after oscillation, whose stability history can be considered as a quasi-periodic response. The Poincare diagram is a number of discrete points, which are concentrated to turnaround in a trajectory area when the speed was 300 r min⁻¹.

4.5 Influence of excitation frequency on linear dynamic characteristics of the system

Taking the 2nd face gear pair as the object, the friction coefficient was 0.01, the input torque was 100 N m⁻¹, the loading speed was 300 r min⁻¹, while the other parameters remained unchanged. The dimensionless excitation frequencies were extracted with 0.1648, 0.0699, 0.0537, 0.0492, 0.0453, and 0.0438 Hz, respectively. The corresponding first six order natural frequencies were 3239.9, 7644.8, 9943.5, 10 852, 11 800, and 12 193 Hz. The displacement curves, phase diagram, and Poincare diagram of the linear dynamic with dimensionless excitation frequencies 0.1648 is shown in Fig. 13 (limited to space, others are not listed).

Figure 13 Displacement, phase and Poincare diagrams of linear system with excitation frequency 0.1648 Hz.

As per Fig. 13, the dimensionless excitation frequency is related to the angular velocity frequency. Similarly, for different frequencies' diagrams, under the condition of inputting a certain speed, as the excitation frequency decreased, the vibration attenuation period of the linear system became longer. The vibration period also became longer after meshing stability, but the vibration displacement amplitude did not change significantly. The phase diagram present stable multiple spiral curves after oscillation, while the Poincare diagram is a number of irregular discrete points.

5 Analysis of nonlinear dynamic characteristics of face gear

The linear model can well approximate the real system dynamics behavior in normal conditions, but overall, it is a simplified model, wherein, it is difficult to obtain a stable and reliable linear approximation (Atanasovska, 2017; Zajicek and Dupal, 2017). Especially, in the analysis of long time-history, the extremely weak nonlinear factors that simplify or even ignore often lead to some unpredictable errors. In general, the simple harmonic motion of the nonlinear system still contains multiple frequency vibration phenomenon with periodicity (or multiple periodicity), but it does not necessarily have the same frequency with the simple harmonic excitation. Under certain excitations, there are several kinds of steady-state vibration characteristics, which make the prediction of the nonlinear dynamic response extremely difficult. Therefore, it is necessary to study the nonlinear dynamics of face gear transmission system, to explore the effect of the main nonlinear factors, in order to reveal the complete dynamic behavior as much as possible. For the multistage face gears transmission system, the teeth surface morphology was reconstructed (or modified) with double crown structure, there were inevitable machining and assembly errors that resulted in a randomly changing gear backlash during the meshing process. Moreover, the double crown tooth surface was difficult for machining, and the frictional response of each tooth surface may cause a certain degree of aperiodic mutation because of the jumping and wear of teeth in the meshing contact area. And it is coupled with time-varying meshing stiffness, damping and external excitation to form a gear dynamics with strong nonlinear characteristics. In this section, the influence of nonlinear dynamic characteristics would be analyzed through the characterization factors such as variable speed stages, friction coefficient, transmission error, excitation frequency, and time-varying stiffness.

With the factors such as the time-variation of the meshing stiffness, the elastic deformation characteristics of the dynamic points contact, the gear backlash, the damping force, the alternating excitation, and the frictional characteristics of the tooth surface, the dimensionless coupled vibration equations can be defined as nonlinear differential equations. Assuming that the dimensionless parameters of the vibration differential equation are: $A_{\mathrm{c}x}=A_{\mathrm{c}z}=\mathrm{0.0261}$, A_dc=0.093, $B_{\mathrm{c}x}=B_{\mathrm{c}z}=\mathrm{0.0809}$, $B_{\mathrm{f}x}=B_{\mathrm{f}z}=\mathrm{0.0071}$, B_dc=0.5509, $A_{\mathrm{f}x}=A_{\mathrm{f}z}=\mathrm{3.2383}\times {\mathrm{10}}^{-\mathrm{4}}$, $A_{\mathrm{df}}=\mathrm{3.8986}\times {\mathrm{10}}^{-\mathrm{4}}$, B_df=0.0231, B_de=0.8806, A_de=0.0149, $C=\mathrm{3.3982}\times {\mathrm{10}}^{-\mathrm{5}}$, and $D=\mathrm{5.1342}\times {\mathrm{10}}^{\mathrm{6}}\times \sin \left(w_i\mathit{\tau }\right)$. Using the 4–5-order Runge–Kutta method, with the friction coefficient as 0.02 and the load as 50 N m⁻¹, the nonlinear characteristics of face gears were analyzed, as shown in Fig. 14.

Figure 14 Displacement, phase, and Poincare diagrams of nonlinear system.

From Fig. 14, under certain load conditions, the response of system is a simple harmonic motion with oscillations in the initial period of the time-history. The vibration displacement of the cylindrical gear was large at startup, and tended to be stable after a period of time. The face gear also presented a relatively stable period simple harmonic. The phase diagram is a plurality of intertwined curves but not crossed, not repeated, and not closed, which is expressed as a quasi-period response. The Poincare diagram is a number of discrete points.

5.1 Influence of variable stage numbers on the nonlinear dynamic characteristics of the system

Other parameters remained unchanged, the teeth number of the face gear was switched to the 1st, 2nd, or 3rd stage, respectively, while the corresponding equivalent mass, rotational inertia, and gear ratio were changed accordingly. The displacement curves, phase diagrams, and Poincare diagrams of the nonlinear dynamics are shown in Figs. 15, 16 and 17.

Figure 15 Displacement, phase, and Poincare diagrams of 1st stage face gear pair nonlinear system.

Figure 16 Displacement, phase, and Poincare diagrams of 2nd stage face gear pair nonlinear system.

Figure 17 Displacement, phase, and Poincare diagrams of 3rd stage face gear pair nonlinear system.

As per Figs. 15, 16, and 17, under the conditions without external excitation, the vibration displacement period of the face gears in all the stages decreased with the increase in the transmission ratio. Because of the gear backlash, the vibration displacement amplitudes of the cylindrical gears generated a certain degree of jumping. The phase diagrams of the nonlinear system is a plurality of intertwined curves but not crossed, not repeated, and not closed, which indicates that there are relatively stable quasi-period characteristics.

5.2 Influence of friction coefficient on the nonlinear dynamic characteristics of the system

Taking 2nd face gear pair as the object, other parameters remained unchanged, the friction coefficient was set to 0.5, 1 and 2 respectively. Because of free idling of the other two stages face gears, but there were still friction damping. Hence, it was necessary to define additional friction coefficients in 1st stage and 3rd stage face gears to be 0.01, and the 2nd stage to be 0.02. The displacement curve, phase diagram, and Poincare diagram of the nonlinear dynamics with friction coefficient 0.5 are shown in Fig. 18 (limited to space, others are not listed).

Figure 18 Displacement, phase, and Poincare diagrams of nonlinear system with friction coefficient 0.5.

From Fig. 18, in the nonlinear system, the friction between the two meshing teeth can induce self-excited vibrations under certain conditions, and because of the gear backlash and frictional damping, the gears generated attenuation free vibration accompanied with contact impact. Similarly, as the friction coefficient became larger, the vibration displacement increasingly deviates from the meshing trajectory. The trend is basically consistent with the linear systems, when the friction coefficient was less than 0.1, the vibration characteristics changed from single cycle to double cycle, and then to quasi-periodic cycle. When greater than 0.1, the vibration displacements in the x and z directions were gradually towards the same direction under the influence of friction damping. The phase diagram is a non-closed curve with multiple spirals. The Poincare diagram is presented as an irregular set of points.

5.3 Influence of transmission error on nonlinear dynamic characteristics of the system

Taking 2nd face gear pair as the object, the friction coefficient was 0.01, the input torque was 100 N m⁻¹, the loading speed was 300 r min⁻¹, while the other parameters remained unchanged. The transmission error fluctuation values were extracted as 0, 0.1, 0.25, and 0.5, respectively. The displacement curve, phase diagram, and Poincare diagram of the nonlinear dynamics with transmission error 0.1 is shown in Fig. 19 (limited to space, others are not listed).

Figure 19 Displacement, phase, and Poincare diagrams of nonlinear system with transmission error fluctuation value 0.1.

As per Fig. 19, when the transmission error was 0.1, the nonlinear system has significant periodic attenuation harmonic motion. The phase diagram is a non-closed spiral curve, and the Poincare diagram is a number of discrete points set concentrating in four regions. Similarly, as the transmission error increased, the vibration displacement generated different degrees of jumping during the simple harmonic motion, whose amplitude also increases accordingly. The phase diagram is a plurality of intertwined and crossed curves combination but neither repeated nor closed. The Poincare diagram mush appears as a number of irregular discrete points.

5.4 Influence of excitation frequency on nonlinear dynamic characteristics of the system

Taking 2nd face gear pair as the object, the friction coefficient was 0.01, the input torque was 100 N m⁻¹, the transmission error fluctuation value was 0.18, while the other parameters remained unchanged. The dimensionless excitation frequency was extracted as 0.1648, 0.2747, 0.5495, and 1.6484 Hz, respectively. The displacement curves, phase diagrams, and Poincare diagrams of the nonlinear dynamics are shown in Figs. 20, 21, 22 and 23.

Figure 20 Displacement, phase, and Poincare diagrams of nonlinear system with excitation frequency 0.1648.

Figure 21 Displacement, phase, and Poincare diagrams of nonlinear system with excitation frequency 0.2747.

Figure 22 Displacement, phase, and Poincare diagrams of nonlinear system with excitation frequency 0.5495.

Figure 23 Displacement, phase, and Poincare diagrams of nonlinear system with excitation frequency 1.6484.

As per Figs. 20, 21, 22 and 23, when the excitation frequency was small, for instance 0.1648, the vibration displacement was a periodic non-harmonic response, which is accompanied with the oscillation attenuation process in the initial period of time-history. When the excitation frequency was 0.2747, the vibration displacement of cylindrical gear reached a stable periodic harmonic response, but face gear still presented a periodic non-harmonic response. The phase diagram is a plurality of non-closed elliptic curves that can be expressed as a quasi-periodic vibration. The Poincare diagram forms a circular area with distinct identification. When the excitation frequency was 0.5495, the vibration displacement of the cylindrical gear presented periodic fluctuation characteristic after oscillation attenuation, and the face gear also tended to harmonic vibration. The phase diagram is a heart-shaped curve region formed by a plurality of non-circular and non-closed curves, which can be expressed as a quasi-periodic vibration response characteristic. The Poincare diagram is a number of discrete points that gradually diverged. When the excitation frequency was 1.6484, the vibration displacements of cylindrical gear and face gear were both the periodic harmonic response, and the phase diagram is a non-circular closed curves set similarly in shape as a curl. The Poincare diagram is a set of irregular discrete points.

6 Conclusion

In this study, single/multistage coupled vibration models were established to derive the coupled vibrations equations. The phonetic relationships of the key parameters such as equivalent mass, rotational inertia, equivalent mesh stiffness, support stiffness, and meshing damping coefficient to dynamic characteristics were studied. The influence rules of the factors such as integrated transmission error, dynamic load, tooth surface friction, loading speed, and load on the transmission system were analyzed. Thus, the linear and nonlinear dynamic characteristics of the coupled vibration differential equations were revealed.

1. The single/multistage coupled vibration models were established using the parameters concentrated method. 7 degrees of freedom vibration equations such as bending vibration, torsional vibration, torsional vibration, and axial vibration were established according to the meshing relationship. The primary characterization parameters such as equivalent mass, integrated meshing stiffness, support stiffness and equivalent damping were defined and described. At the same time, the conditions and the main inducing factors of internal excitation or external excitation of the transmission system were also clarified. Thus, a dynamic model suitable for the multistage face gears transmission has been established and optimized.

2. The vibration equation was processed by reducing the degrees of freedom, and the corresponding parameters were subjected to dimensionless normalization. The ODE45 integral method was introduced to solve the vibration equations. The linear and nonlinear vibration conditions were defined, and the effects of various parameters or excitations of the linear or nonlinear system were analyzed in detail to accurately solve the vibration performance and power transmission performance of the transmission system.

3. The linear dynamic characteristics can approximate the vibration response of the real system, but it can only be expressed as a simplified model. It is difficult to maintain the stability and reliability under the influence of various nonlinear errors, and its vibration characteristics are also difficult to predict. Therefore, the accuracy of the analysis results for the dynamic characteristics of the transmission system can be improved by comprehensively synthesizing the interaction of nonlinear factors such as external excitation and internal excitation.