Journal topic
Mech. Sci., 10, 475–495, 2019
https://doi.org/10.5194/ms-10-475-2019
Mech. Sci., 10, 475–495, 2019
https://doi.org/10.5194/ms-10-475-2019

Research article 24 Sep 2019

Research article | 24 Sep 2019

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

Numerical modeling and dynamic characteristics study of coupling vibration of multistage face gearsplanetary transmission
Xingbin Chen1,2,3, Qingchun Hu1, Zhongyang Xu1, and Chune Zhu2,4 Xingbin Chen et al.
• 1School of Mechanical and Automotive engineering,South China University of Technology, Guangzhou, 510641, China
• 2Guangzhou Mechanical Engineering Research Institute Co., Ltd, Guangzhou, 510700, China
• 3CRAT Testing & Certification Company Ltd., Guangzhou, 510700, China
• 4Institute for Biomedical and Pharmaceutical Sciences, Guangdong University of Technology, Guangzhou, 510006, China

Correspondence: Qingchun Hu (huqc@scut.edu.cn)

Abstract

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

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. Kcx and Kcz are the stiffness coefficients of the cylindrical gear in the x and z directions, respectively. Kfx and Kfy are the stiffness coefficients of the face gear in the x and z directions, respectively. Kh is the stiffness coefficient of the gear pair. Ccx and Ccy are the damping coefficients of the cylindrical gear in the x and z directions, respectively. Cfx and Cfy are the damping coefficients of the face gear in the x and y directions, respectively. Ch is the damping coefficient of the gear pair. eh is the static integrated transmission error of the gear pair. bh 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. Tc and Tf are the torques of cylindrical gear and face gear, respectively. Cf12 and Cf23 are the damping coefficients of friction and backlash between the stages of multistage face gears.

Figure 1Kinematics model of single-stage face gear.

Figure 2Kinematics 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 xc and zc. The translation of face gear in x axis is xf.

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 zf.

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}\stackrel{\mathrm{¨}}{q}+\mathbf{C}\stackrel{\mathrm{˙}}{q}+\mathbf{K}q+{\mathit{\varsigma }}_{q}^{\mathrm{T}}\mathit{\lambda }=\mathbf{F}+\mathbf{Q}\\ \mathit{\varsigma }\left(q,t\right)=\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 Fn, and the component forces Fx and Fz 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}}{\stackrel{\mathrm{˙}}{d}}_{\mathrm{n}}\\ {\mathbf{F}}_{x}={\mathbf{F}}_{\mathrm{n}}\mathrm{cos}{\mathit{\alpha }}_{\mathrm{n}}\\ {\mathbf{F}}_{z}={\mathbf{F}}_{\mathrm{n}}\mathrm{sin}{\mathit{\alpha }}_{\mathrm{n}}\end{array}\right\\end{array}$

In Eq. (3), km(t) is the time-varying meshing stiffness, f(dn) is the backlash function, cm is the meshing damping, and dn 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}}& =\left({x}_{\mathrm{c}}+{r}_{\mathrm{cm}}{\mathit{\theta }}_{\mathrm{c}y}\right)\mathrm{cos}{\mathit{\alpha }}_{\mathrm{n}}+{z}_{\mathrm{c}}\mathrm{sin}{\mathit{\alpha }}_{\mathrm{n}}-\left({x}_{\mathrm{f}}+{r}_{\mathrm{fm}}{\mathit{\theta }}_{\mathrm{f}z}\right)\\ & \mathrm{cos}{\mathit{\alpha }}_{\mathrm{n}}-\left({z}_{\mathrm{f}}+{r}_{\mathrm{fm}}{\mathit{\theta }}_{\mathrm{f}x}\right)\mathrm{sin}{\mathit{\alpha }}_{\mathrm{n}}-{e}_{\mathrm{n}}\left(t\right)\end{array}\end{array}$

In Eq. (4), rcm and rfm are the meshing point radii of the cylindrical gear and face gear, respectively. αn are the normal pressure angles of meshing points. en(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}}{\stackrel{\mathrm{¨}}{x}}_{\mathrm{c}}+{c}_{\mathrm{c}x}{\stackrel{\mathrm{˙}}{x}}_{\mathrm{c}}+{k}_{\mathrm{c}x}{x}_{\mathrm{c}}=-{\mathbf{F}}_{x}-\mathit{\mu }{\overline{f}}_{\mathit{\mu }}{\mathbf{F}}_{\mathrm{n}}\mathrm{sin}{\mathit{\alpha }}_{\mathrm{n}}\\ {\mathbf{M}}_{\mathrm{c}}{\stackrel{\mathrm{¨}}{z}}_{\mathrm{c}}+{c}_{\mathrm{c}z}{\stackrel{\mathrm{˙}}{z}}_{\mathrm{c}}+{k}_{\mathrm{c}z}{z}_{\mathrm{c}}=-{\mathbf{F}}_{z}-\mathit{\mu }{\overline{f}}_{\mathit{\mu }}{\mathbf{F}}_{\mathrm{n}}\mathrm{cos}{\mathit{\alpha }}_{\mathrm{n}}\\ {\mathbf{J}}_{\mathrm{c}y}{\stackrel{\mathrm{¨}}{\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}}\mathrm{sin}{\mathit{\alpha }}_{\mathrm{n}}\\ {\mathbf{M}}_{\mathrm{f}}{\stackrel{\mathrm{¨}}{x}}_{\mathrm{f}}+{c}_{\mathrm{f}x}{\stackrel{\mathrm{˙}}{x}}_{\mathrm{f}}+{k}_{\mathrm{f}x}{x}_{\mathrm{f}}={\mathbf{F}}_{x}+\mathit{\mu }{\overline{f}}_{\mathit{\mu }}{\mathbf{F}}_{\mathrm{n}}\mathrm{sin}{\mathit{\alpha }}_{\mathrm{n}}\\ {\mathbf{M}}_{\mathrm{f}}{\stackrel{\mathrm{¨}}{z}}_{\mathrm{f}}+{c}_{\mathrm{f}z}{\stackrel{\mathrm{˙}}{z}}_{\mathrm{f}}+{k}_{\mathrm{f}z}{z}_{\mathrm{f}}={\mathbf{F}}_{z}+\mathit{\mu }{\overline{f}}_{\mathit{\mu }}{\mathbf{F}}_{\mathrm{n}}\mathrm{cos}{\mathit{\alpha }}_{\mathrm{n}}\\ {\mathbf{J}}_{\mathrm{f}x}{\stackrel{\mathrm{¨}}{\mathit{\theta }}}_{\mathrm{f}x}+{c}_{{\mathit{\theta }}_{\mathrm{f}x}}{\stackrel{\mathrm{˙}}{\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}}\mathrm{sin}{\mathit{\alpha }}_{\mathrm{n}}\\ {\mathbf{J}}_{\mathrm{f}z}{\stackrel{\mathrm{¨}}{\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}}\mathrm{cos}{\mathit{\alpha }}_{\mathrm{n}}\end{array}\right\\end{array}$

In Eq. (5), Mc and Mf are the concentrated mass of cylindrical gear and face gear, respectively. kij and ${c}_{ij}\left(i=c,f;j=x,z\right)$ are the support stiffness and damping of the cylindrical gear or face gear along the x and z axes, respectively. Jcy is the rotational inertia of the cylindrical gear along the y axes. Jfx and Jfz are the rotational inertias of face gear in the x and z axes, respectively. Tc is the input torque acting on the cylindrical gear. Tf 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. rcb and rca are the base radius and addendum radius of the cylindrical gear, respectively. ε is the overlap ratio of the gear pair. scb is the base tooth pitch. nc is the speed of the cylindrical gear. ${\overline{f}}_{\mathit{\mu }}=\text{sign}\left({L}_{A}-{r}_{\mathrm{cb}}\mathrm{tan}{\mathit{\alpha }}_{\mathrm{n}}\right)$ 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}{\stackrel{\mathrm{¨}}{\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}}\mathrm{sin}{\mathit{\alpha }}_{\mathrm{n}}\\ {\mathbf{J}}_{\mathrm{f}z}{\stackrel{\mathrm{¨}}{\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}}\mathrm{cos}{\mathit{\alpha }}_{\mathrm{n}}\end{array}\right\\end{array}$

After rewriting:

$\begin{array}{}\text{(7)}& \left\{\begin{array}{l}{r}_{\mathrm{cm}}{\stackrel{\mathrm{¨}}{\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}}\mathrm{sin}{\mathit{\alpha }}_{\mathrm{n}}}{{\mathbf{J}}_{\mathrm{c}y}}{r}_{\mathrm{cm}}\\ {r}_{\mathrm{fm}}{\stackrel{\mathrm{¨}}{\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}}\mathrm{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}}{\stackrel{\mathrm{¨}}{\mathit{\theta }}}_{\mathrm{c}y}-{r}_{\mathrm{fm}}{\stackrel{\mathrm{¨}}{\mathit{\theta }}}_{\mathrm{f}z}& =\frac{{\mathbf{T}}_{\mathrm{c}}-{\mathbf{F}}_{\mathrm{n}}{r}_{\mathrm{cm}}}{{\mathbf{J}}_{\mathrm{c}y}}{r}_{\mathrm{cm}}-\frac{{\mathbf{F}}_{\mathrm{n}}{r}_{\mathrm{fm}}-{\mathbf{T}}_{\mathrm{f}}}{{\mathbf{J}}_{\mathrm{f}z}}{r}_{\mathrm{fm}}\\ & +\frac{\mathit{\mu }{\overline{f}}_{\mathit{\mu }}{\mathbf{L}}_{A}{\mathbf{F}}_{\mathrm{n}}\mathrm{sin}{\mathit{\alpha }}_{\mathrm{n}}}{{\mathbf{J}}_{\mathrm{c}y}}{r}_{\mathrm{cm}}\\ & -\frac{\mathit{\mu }{\overline{f}}_{\mathit{\mu }}{r}_{\mathrm{fm}}{\mathbf{F}}_{\mathrm{n}}\mathrm{cos}{\mathit{\alpha }}_{\mathrm{n}}}{{\mathbf{J}}_{\mathrm{f}z}}{r}_{\mathrm{fm}}\end{array}\end{array}$

The first and second derivative ${\stackrel{\mathrm{˙}}{d}}_{\mathrm{n}}$, ${\stackrel{\mathrm{¨}}{d}}_{\mathrm{n}}$ of displacement dn in the normal direction of mesh points are introduced:

$\begin{array}{}\text{(9)}& \left\{\begin{array}{ll}{\stackrel{\mathrm{˙}}{d}}_{\mathrm{n}}& =\left({\stackrel{\mathrm{˙}}{x}}_{\mathrm{c}}-{\stackrel{\mathrm{˙}}{x}}_{\mathrm{f}}\right)\mathrm{cos}{\mathit{\alpha }}_{\mathrm{n}}+\left({r}_{\mathrm{cm}}{\stackrel{\mathrm{˙}}{\mathit{\theta }}}_{\mathrm{c}y}-{r}_{\mathrm{fm}}{\stackrel{\mathrm{˙}}{\mathit{\theta }}}_{\mathrm{f}z}\right)\mathrm{cos}{\mathit{\alpha }}_{\mathrm{n}}\\ & +\left({\stackrel{\mathrm{˙}}{z}}_{\mathrm{c}}-{\stackrel{\mathrm{˙}}{z}}_{\mathrm{f}}-{r}_{\mathrm{fm}}{\stackrel{\mathrm{˙}}{\mathit{\theta }}}_{\mathrm{f}x}\right)\mathrm{sin}{\mathit{\alpha }}_{\mathrm{n}}-{\stackrel{\mathrm{˙}}{e}}_{\mathrm{n}}\left(t\right)\\ {\stackrel{\mathrm{¨}}{d}}_{\mathrm{n}}& =\left({\stackrel{\mathrm{¨}}{x}}_{\mathrm{c}}-{\stackrel{\mathrm{¨}}{x}}_{\mathrm{f}}\right)\mathrm{cos}{\mathit{\alpha }}_{\mathrm{n}}+\left({r}_{\mathrm{cm}}{\stackrel{\mathrm{¨}}{\mathit{\theta }}}_{\mathrm{c}y}-{r}_{\mathrm{fm}}{\stackrel{\mathrm{¨}}{\mathit{\theta }}}_{\mathrm{f}z}\right)\mathrm{cos}{\mathit{\alpha }}_{\mathrm{n}}\\ & +\left({\stackrel{\mathrm{¨}}{z}}_{\mathrm{c}}-{\stackrel{\mathrm{¨}}{z}}_{\mathrm{f}}-{r}_{\mathrm{fm}}{\stackrel{\mathrm{¨}}{\mathit{\theta }}}_{\mathrm{f}x}\right)\mathrm{sin}{\mathit{\alpha }}_{\mathrm{n}}-{\stackrel{\mathrm{¨}}{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}{\stackrel{\mathrm{¨}}{d}}_{\mathrm{n}}& =\left({\stackrel{\mathrm{¨}}{x}}_{\mathrm{c}}-{\stackrel{\mathrm{¨}}{x}}_{\mathrm{f}}\right)\mathrm{cos}{\mathit{\alpha }}_{\mathrm{n}}+\left(\frac{{\mathbf{T}}_{\mathrm{c}}-{\mathbf{F}}_{\mathrm{n}}{r}_{\mathrm{cm}}}{{\mathbf{J}}_{\mathrm{c}y}}{r}_{\mathrm{cm}}-\frac{{\mathbf{F}}_{\mathrm{n}}{r}_{\mathrm{fm}}-{\mathbf{T}}_{\mathrm{f}}}{{\mathbf{J}}_{\mathrm{f}z}}{r}_{\mathrm{fm}}\right\\ & +\frac{\mathit{\mu }{\overline{f}}_{\mathit{\mu }}{\mathbf{L}}_{A}{\mathbf{F}}_{\mathrm{n}}\mathrm{sin}{\mathit{\alpha }}_{\mathrm{n}}}{{\mathbf{J}}_{\mathrm{c}y}}{r}_{\mathrm{cm}}-\frac{\mathit{\mu }{\overline{f}}_{\mathit{\mu }}{r}_{\mathrm{fm}}{\mathbf{F}}_{\mathrm{n}}\mathrm{cos}{\mathit{\alpha }}_{\mathrm{n}}}{{\mathbf{J}}_{\mathrm{f}z}}{r}_{\mathrm{fm}})\mathrm{cos}{\mathit{\alpha }}_{\mathrm{n}}\\ & +\left({\stackrel{\mathrm{¨}}{z}}_{\mathrm{c}}-{\stackrel{\mathrm{¨}}{z}}_{\mathrm{f}}-{r}_{\mathrm{fm}}{\stackrel{\mathrm{¨}}{\mathit{\theta }}}_{\mathrm{f}x}\right)\mathrm{sin}{\mathit{\alpha }}_{\mathrm{n}}-{\stackrel{\mathrm{¨}}{e}}_{\mathrm{n}}\left(t\right)\end{array}\end{array}$

Me is substituted and rearranged as:

$\begin{array}{}\text{(11)}& \begin{array}{rl}{\mathbf{M}}_{e}{\stackrel{\mathrm{¨}}{d}}_{\mathrm{n}}& -{\mathbf{M}}_{e}\left({\stackrel{\mathrm{¨}}{x}}_{\mathrm{c}}-{\stackrel{\mathrm{¨}}{x}}_{\mathrm{f}}\right)\mathrm{cos}{\mathit{\alpha }}_{\mathrm{n}}-{\mathbf{M}}_{e}\left({\stackrel{\mathrm{¨}}{z}}_{\mathrm{c}}-{\stackrel{\mathrm{¨}}{z}}_{\mathrm{f}}\right)\mathrm{sin}{\mathit{\alpha }}_{\mathrm{n}}\\ & +{\mathbf{M}}_{e}{\stackrel{\mathrm{¨}}{e}}_{\mathrm{n}}\left(t\right)=\frac{{\mathbf{M}}_{e}{\mathbf{T}}_{\mathrm{c}}{r}_{\mathrm{cm}}\mathrm{cos}{\mathit{\alpha }}_{\mathrm{n}}}{{\mathbf{J}}_{\mathrm{c}y}}+\frac{{\mathbf{M}}_{e}{\mathbf{T}}_{\mathrm{f}}{r}_{\mathrm{fm}}\mathrm{cos}{\mathit{\alpha }}_{\mathrm{n}}}{{\mathbf{J}}_{\mathrm{f}z}}\\ & +\frac{{\mathbf{M}}_{e}\mathit{\mu }{\overline{f}}_{\mathit{\mu }}{\mathbf{L}}_{A}{r}_{\mathrm{cm}}{\mathbf{F}}_{\mathrm{n}}\mathrm{sin}{\mathit{\alpha }}_{\mathrm{n}}\mathrm{cos}{\mathit{\alpha }}_{\mathrm{n}}}{{\mathbf{J}}_{\mathrm{c}y}}\\ & -\frac{{\mathbf{M}}_{e}\mathit{\mu }{\overline{f}}_{\mathit{\mu }}{r}_{\mathrm{fm}}^{\mathrm{2}}{\mathbf{F}}_{\mathrm{n}}{\mathrm{cos}}^{\mathrm{2}}{\mathit{\alpha }}_{\mathrm{n}}}{{\mathbf{J}}_{\mathrm{f}z}}\\ & -\mathrm{cos}{\mathit{\alpha }}_{\mathrm{n}}{\mathbf{F}}_{\mathrm{n}}-{\mathbf{M}}_{e}{r}_{\mathrm{fm}}{\stackrel{\mathrm{¨}}{\mathit{\theta }}}_{\mathrm{f}x}\mathrm{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}}\mathrm{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}}\mathrm{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}}\mathrm{sin}{\mathit{\alpha }}_{\mathrm{n}}\mathrm{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}}{\mathrm{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}{\stackrel{\mathrm{¨}}{d}}_{\mathrm{n}}& -{\mathbf{M}}_{e}\left({\stackrel{\mathrm{¨}}{x}}_{\mathrm{c}}-{\stackrel{\mathrm{¨}}{x}}_{\mathrm{f}}\right)\mathrm{cos}{\mathit{\alpha }}_{\mathrm{n}}-{\mathbf{M}}_{e}\left({\stackrel{\mathrm{¨}}{z}}_{\mathrm{c}}-{\stackrel{\mathrm{¨}}{z}}_{\mathrm{f}}\right)\mathrm{sin}{\mathit{\alpha }}_{\mathrm{n}}+{\mathbf{M}}_{e}{\stackrel{\mathrm{¨}}{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 }-\mathrm{cos}{\mathit{\alpha }}_{\mathrm{n}}{\mathbf{F}}_{\mathrm{n}}-{\mathbf{M}}_{e}{r}_{\mathrm{fm}}{\stackrel{\mathrm{¨}}{\mathit{\theta }}}_{\mathrm{f}x}\mathrm{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}}{\stackrel{\mathrm{¨}}{x}}_{\mathrm{c}}+{c}_{\mathrm{c}x}{\stackrel{\mathrm{˙}}{x}}_{\mathrm{c}}+{k}_{\mathrm{c}x}{x}_{\mathrm{c}}=-{\mathbf{F}}_{x}-\mathit{\mu }{\overline{f}}_{\mathit{\mu }}{\mathbf{F}}_{\mathrm{n}}\mathrm{sin}{\mathit{\alpha }}_{\mathrm{n}}\\ {\mathbf{M}}_{\mathrm{c}}{\stackrel{\mathrm{¨}}{z}}_{\mathrm{c}}+{c}_{\mathrm{c}z}{\stackrel{\mathrm{˙}}{z}}_{\mathrm{c}}+{k}_{\mathrm{c}z}{z}_{\mathrm{c}}=-{\mathbf{F}}_{z}-\mathit{\mu }{\overline{f}}_{\mathit{\mu }}{\mathbf{F}}_{\mathrm{n}}\mathrm{cos}{\mathit{\alpha }}_{\mathrm{n}}\\ {\mathbf{M}}_{\mathrm{f}}{\stackrel{\mathrm{¨}}{x}}_{\mathrm{f}}+{c}_{\mathrm{f}x}{\stackrel{\mathrm{˙}}{x}}_{\mathrm{f}}+{k}_{\mathrm{f}x}{x}_{\mathrm{f}}={\mathbf{F}}_{x}+\mathit{\mu }{\overline{f}}_{\mathit{\mu }}{\mathbf{F}}_{\mathrm{n}}\mathrm{sin}{\mathit{\alpha }}_{\mathrm{n}}\\ {\mathbf{M}}_{\mathrm{f}}{\stackrel{\mathrm{¨}}{z}}_{\mathrm{f}}+{c}_{\mathrm{f}z}{\stackrel{\mathrm{˙}}{z}}_{\mathrm{f}}+{k}_{\mathrm{f}z}{z}_{\mathrm{f}}={\mathbf{F}}_{z}+\mathit{\mu }{\overline{f}}_{\mathit{\mu }}{\mathbf{F}}_{\mathrm{n}}\mathrm{cos}{\mathit{\alpha }}_{\mathrm{n}}\\ {\mathbf{J}}_{\mathrm{f}x}{\stackrel{\mathrm{¨}}{\mathit{\theta }}}_{\mathrm{f}x}+{c}_{{\mathit{\theta }}_{\mathrm{f}x}}{\stackrel{\mathrm{˙}}{\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}}\mathrm{sin}{\mathit{\alpha }}_{\mathrm{n}}\\ {\mathbf{M}}_{e}{\stackrel{\mathrm{¨}}{d}}_{\mathrm{n}}-{\mathbf{M}}_{e}\left({\stackrel{\mathrm{¨}}{x}}_{\mathrm{c}}-{\stackrel{\mathrm{¨}}{x}}_{\mathrm{f}}\right)\mathrm{cos}{\mathit{\alpha }}_{\mathrm{n}}-{\mathbf{M}}_{e}\left({\stackrel{\mathrm{¨}}{z}}_{\mathrm{c}}-{\stackrel{\mathrm{¨}}{z}}_{\mathrm{f}}\right)\mathrm{sin}{\mathit{\alpha }}_{\mathrm{n}}+{\mathbf{M}}_{e}{\stackrel{\mathrm{¨}}{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 }-\mathrm{cos}{\mathit{\alpha }}_{\mathrm{n}}{\mathbf{F}}_{\mathrm{n}}-{\mathbf{M}}_{e}{r}_{\mathrm{fm}}{\stackrel{\mathrm{¨}}{\mathit{\theta }}}_{\mathrm{f}x}\mathrm{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}}{\stackrel{\mathrm{¨}}{x}}_{\mathrm{c}}+{c}_{\mathrm{c}x}{\stackrel{\mathrm{˙}}{x}}_{\mathrm{c}}+{k}_{\mathrm{c}x}{x}_{\mathrm{c}}=-\left(\mathrm{cos}{\mathit{\alpha }}_{\mathrm{n}}+\mathit{\mu }{\overline{f}}_{\mathit{\mu }}\mathrm{sin}{\mathit{\alpha }}_{\mathrm{n}}\right){\mathbf{F}}_{\mathrm{n}}\\ {\mathbf{M}}_{\mathrm{c}}{\stackrel{\mathrm{¨}}{z}}_{\mathrm{c}}+{c}_{\mathrm{c}z}{\stackrel{\mathrm{˙}}{z}}_{\mathrm{c}}+{k}_{\mathrm{c}z}{z}_{\mathrm{c}}=-\left(\mathrm{sin}{\mathit{\alpha }}_{\mathrm{n}}-\mathit{\mu }{\overline{f}}_{\mathit{\mu }}\mathrm{cos}{\mathit{\alpha }}_{\mathrm{n}}\right){\mathbf{F}}_{\mathrm{n}}\\ {\mathbf{M}}_{\mathrm{f}}{\stackrel{\mathrm{¨}}{x}}_{\mathrm{f}}+{c}_{\mathrm{f}x}{\stackrel{\mathrm{˙}}{x}}_{\mathrm{f}}+{k}_{\mathrm{f}x}{x}_{\mathrm{f}}=\left(\mathrm{cos}{\mathit{\alpha }}_{\mathrm{n}}+\mathit{\mu }{\overline{f}}_{\mathit{\mu }}\mathrm{sin}{\mathit{\alpha }}_{\mathrm{n}}\right){\mathbf{F}}_{\mathrm{n}}\\ {\mathbf{M}}_{\mathrm{f}}{\stackrel{\mathrm{¨}}{z}}_{\mathrm{f}}+{c}_{\mathrm{f}z}{\stackrel{\mathrm{˙}}{z}}_{\mathrm{f}}+{k}_{\mathrm{f}z}{z}_{\mathrm{f}}=\left(\mathrm{sin}{\mathit{\alpha }}_{\mathrm{n}}+\mathit{\mu }{\overline{f}}_{\mathit{\mu }}\mathrm{cos}{\mathit{\alpha }}_{\mathrm{n}}\right){\mathbf{F}}_{\mathrm{n}}\\ {\mathbf{M}}_{e}{\stackrel{\mathrm{¨}}{d}}_{\mathrm{n}}-{\mathbf{M}}_{e}\left({\stackrel{\mathrm{¨}}{x}}_{\mathrm{c}}-{\stackrel{\mathrm{¨}}{x}}_{\mathrm{f}}\right)\mathrm{cos}{\mathit{\alpha }}_{\mathrm{n}}-{\mathbf{M}}_{e}\left({\stackrel{\mathrm{¨}}{z}}_{\mathrm{c}}-{\stackrel{\mathrm{¨}}{z}}_{\mathrm{f}}\right)\mathrm{sin}{\mathit{\alpha }}_{\mathrm{n}}\\ +{\mathbf{M}}_{e}{\stackrel{\mathrm{¨}}{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 }-\mathrm{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 ND 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}}\\ & -\left(\mathrm{cos}{\mathit{\alpha }}_{\mathrm{n}}+\mathit{\mu }{\overline{f}}_{\mathit{\mu }}\mathrm{sin}{\mathit{\alpha }}_{\mathrm{n}}\right)\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)\\ & -\left(\mathrm{cos}{\mathit{\alpha }}_{\mathrm{n}}+\mathit{\mu }{\overline{f}}_{\mathit{\mu }}\mathrm{sin}{\mathit{\alpha }}_{\mathrm{n}}\right)\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}}\\ & +\left(\mathrm{sin}{\mathit{\alpha }}_{\mathrm{n}}-\mathit{\mu }{\overline{f}}_{\mathit{\mu }}\mathrm{cos}{\mathit{\alpha }}_{\mathrm{n}}\right)\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)\\ & +\left(\mathrm{sin}{\mathit{\alpha }}_{\mathrm{n}}-\mathit{\mu }{\overline{f}}_{\mathit{\mu }}\mathrm{cos}{\mathit{\alpha }}_{\mathrm{n}}\right)\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}}\\ & -\left(\mathrm{cos}{\mathit{\alpha }}_{\mathrm{n}}+\mathit{\mu }{\overline{f}}_{\mathit{\mu }}\mathrm{sin}{\mathit{\alpha }}_{\mathrm{n}}\right)\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)\\ & -\left(\mathrm{cos}{\mathit{\alpha }}_{\mathrm{n}}+\mathit{\mu }{\overline{f}}_{\mathit{\mu }}\mathrm{sin}{\mathit{\alpha }}_{\mathrm{n}}\right)\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}}\\ & +\left(\mathrm{sin}{\mathit{\alpha }}_{\mathrm{n}}-\mathit{\mu }{\overline{f}}_{\mathit{\mu }}\mathrm{cos}{\mathit{\alpha }}_{\mathrm{n}}\right)\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)\\ & +\left(\mathrm{sin}{\mathit{\alpha }}_{\mathrm{n}}-\mathit{\mu }{\overline{f}}_{\mathit{\mu }}\mathrm{cos}{\mathit{\alpha }}_{\mathrm{n}}\right)\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)\mathrm{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)\mathrm{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)\mathrm{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)\mathrm{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}}}{\stackrel{\mathrm{¨}}{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}{\stackrel{\mathrm{˙}}{e}}_{\mathrm{nf}}\left(\mathit{\tau }\right)={e}_{r}{\mathit{\omega }}_{i}\mathrm{cos}\left({\mathit{\omega }}_{i}\mathit{\tau }+{\mathit{\phi }}_{ipa}\right)\\ {\stackrel{\mathrm{¨}}{e}}_{\mathrm{n}f}\left(\mathit{\tau }\right)=-{e}_{r}{\mathit{\omega }}_{i}^{\mathrm{2}}\mathrm{sin}\left({\mathit{\omega }}_{i}\mathit{\tau }+{\mathit{\phi }}_{ipa}\right)\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}}\mathrm{sin}\left({\mathit{\omega }}_{i}\mathit{\tau }+{\mathit{\phi }}_{ipa}\right)}{{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}}\\ & -\left(\mathrm{cos}{\mathit{\alpha }}_{\mathrm{n}}+\mathit{\mu }{\overline{f}}_{\mathit{\mu }}\mathrm{sin}{\mathit{\alpha }}_{\mathrm{n}}\right){B}_{\mathrm{dc}}{f}_{\mathrm{b}}\left({d}_{\mathrm{d}}\right)\\ & -\mathrm{2}\left(\mathrm{cos}{\mathit{\alpha }}_{\mathrm{n}}+\mathit{\mu }{\overline{f}}_{\mathit{\mu }}\mathrm{sin}{\mathit{\alpha }}_{\mathrm{n}}\right){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}}\\ & +\left(\mathrm{sin}{\mathit{\alpha }}_{\mathrm{n}}-\mathit{\mu }{\overline{f}}_{\mathit{\mu }}\mathrm{cos}{\mathit{\alpha }}_{\mathrm{n}}\right){B}_{\mathrm{dc}}{f}_{\mathrm{b}}\left({d}_{\mathrm{d}}\right)\\ & +\mathrm{2}\left(\mathrm{sin}{\mathit{\alpha }}_{\mathrm{n}}-\mathit{\mu }{\overline{f}}_{\mathit{\mu }}\mathrm{cos}{\mathit{\alpha }}_{\mathrm{n}}\right){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}}\\ & -\left(\mathrm{cos}{\mathit{\alpha }}_{\mathrm{n}}+\mathit{\mu }{\overline{f}}_{\mathit{\mu }}\mathrm{sin}{\mathit{\alpha }}_{\mathrm{n}}\right){B}_{\mathrm{df}}{f}_{\mathrm{b}}\left({d}_{\mathrm{d}}\right)\\ & -\mathrm{2}\left(\mathrm{cos}{\mathit{\alpha }}_{\mathrm{n}}+\mathit{\mu }{\overline{f}}_{\mathit{\mu }}\mathrm{sin}{\mathit{\alpha }}_{\mathrm{n}}\right){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}}\\ & +\left(\mathrm{sin}{\mathit{\alpha }}_{\mathrm{n}}-\mathit{\mu }{\overline{f}}_{\mathit{\mu }}\mathrm{cos}{\mathit{\alpha }}_{\mathrm{n}}\right){B}_{\mathrm{df}}{f}_{\mathrm{b}}\left({d}_{\mathrm{d}}\right)\\ & +\mathrm{2}\left(\mathrm{sin}{\mathit{\alpha }}_{\mathrm{n}}-\mathit{\mu }{\overline{f}}_{\mathit{\mu }}\mathrm{cos}{\mathit{\alpha }}_{\mathrm{n}}\right){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}}}\mathrm{cos}{\mathit{\alpha }}_{\mathrm{n}}+\frac{{d}^{\mathrm{2}}{x}_{\mathrm{df}}}{\mathrm{d}{\mathit{\tau }}^{\mathrm{2}}}\mathrm{cos}{\mathit{\alpha }}_{\mathrm{n}}-\frac{{d}^{\mathrm{2}}{z}_{\mathrm{dc}}}{\mathrm{d}{\mathit{\tau }}^{\mathrm{2}}}\mathrm{sin}{\mathit{\alpha }}_{\mathrm{n}}\\ & +\frac{{d}^{\mathrm{2}}{z}_{\mathrm{df}}}{\mathrm{d}{\mathit{\tau }}^{\mathrm{2}}}\mathrm{sin}{\mathit{\alpha }}_{\mathrm{n}}+{B}_{\mathrm{d}e}\mathrm{cos}{\mathit{\alpha }}_{\mathrm{n}}{f}_{\mathrm{b}}\left({d}_{\mathrm{d}}\right)\\ & +\mathrm{2}{A}_{\mathrm{d}e}\mathrm{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{\left\{}{\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{\right\}}}^{\mathrm{T}}\\ & =\mathit{\left\{}{x}_{\mathrm{dc}},{\stackrel{\mathrm{˙}}{x}}_{\mathrm{dc}},{z}_{\mathrm{dc}},{\stackrel{\mathrm{˙}}{z}}_{\mathrm{dc}},{x}_{\mathrm{df}},{\stackrel{\mathrm{˙}}{x}}_{\mathrm{df}},{z}_{\mathrm{df}},{\stackrel{\mathrm{˙}}{z}}_{\mathrm{df}},{d}_{\mathrm{d}},{\stackrel{\mathrm{˙}}{d}}_{\mathrm{d}}{\mathit{\right\}}}^{\mathrm{T}}\end{array}\end{array}$

which are:

$\begin{array}{}\text{(21)}& \left\{\begin{array}{ll}{\stackrel{\mathrm{˙}}{\mathit{\lambda }}}_{\mathrm{1}}& ={\mathit{\lambda }}_{\mathrm{2}}\\ {\stackrel{\mathrm{˙}}{\mathit{\lambda }}}_{\mathrm{2}}& =-\mathrm{2}{A}_{\mathrm{c}x}{\mathit{\lambda }}_{\mathrm{2}}-{B}_{\mathrm{c}x}{\mathit{\lambda }}_{\mathrm{1}}+\left(\mathrm{cos}{\mathit{\alpha }}_{\mathrm{n}}\\ & +\mathit{\mu }{\overline{f}}_{\mathit{\mu }}\mathrm{sin}{\mathit{\alpha }}_{\mathrm{n}}\right){B}_{\mathrm{dc}}{f}_{\mathrm{b}}\left({\mathit{\lambda }}_{\mathrm{9}}\right)\\ & +\mathrm{2}\left(\mathrm{cos}{\mathit{\alpha }}_{\mathrm{n}}+\mathit{\mu }{\overline{f}}_{\mathit{\mu }}\mathrm{sin}{\mathit{\alpha }}_{\mathrm{n}}\right){A}_{\mathrm{dc}}{\mathit{\lambda }}_{\mathrm{10}}\\ {\stackrel{\mathrm{˙}}{\mathit{\lambda }}}_{\mathrm{3}}& ={\mathit{\lambda }}_{\mathrm{4}}\\ {\stackrel{\mathrm{˙}}{\mathit{\lambda }}}_{\mathrm{4}}& =-\mathrm{2}{A}_{\mathrm{c}z}{\mathit{\lambda }}_{\mathrm{4}}-{B}_{\mathrm{c}z}{\mathit{\lambda }}_{\mathrm{3}}\\ & -\left(\mathrm{sin}{\mathit{\alpha }}_{\mathrm{n}}-\mathit{\mu }{\overline{f}}_{\mathit{\mu }}\mathrm{cos}{\mathit{\alpha }}_{\mathrm{n}}\right){B}_{\mathrm{dc}}{f}_{\mathrm{b}}\left({\mathit{\lambda }}_{\mathrm{9}}\right)\\ & -\mathrm{2}\left(\mathrm{sin}{\mathit{\alpha }}_{\mathrm{n}}-\mathit{\mu }{\overline{f}}_{\mathit{\mu }}\mathrm{cos}{\mathit{\alpha }}_{\mathrm{n}}\right){A}_{\mathrm{dc}}{\mathit{\lambda }}_{\mathrm{10}}\\ {\stackrel{\mathrm{˙}}{\mathit{\lambda }}}_{\mathrm{5}}& ={\mathit{\lambda }}_{\mathrm{6}}\\ {\stackrel{\mathrm{˙}}{\mathit{\lambda }}}_{\mathrm{6}}& =-\mathrm{2}{A}_{\mathrm{f}x}{\mathit{\lambda }}_{\mathrm{6}}-{B}_{\mathrm{f}x}{\mathit{\lambda }}_{\mathrm{5}}\\ & +\left(\mathrm{cos}{\mathit{\alpha }}_{\mathrm{n}}+\mathit{\mu }{\overline{f}}_{\mathit{\mu }}\mathrm{sin}{\mathit{\alpha }}_{\mathrm{n}}\right){B}_{\mathrm{df}}{f}_{\mathrm{b}}\left({\mathit{\lambda }}_{\mathrm{9}}\right)\\ & +\mathrm{2}\left(\mathrm{cos}{\mathit{\alpha }}_{\mathrm{n}}+\mathit{\mu }{\overline{f}}_{\mathit{\mu }}\mathrm{sin}{\mathit{\alpha }}_{\mathrm{n}}\right){A}_{\mathrm{df}}{\mathit{\lambda }}_{\mathrm{10}}\\ {\stackrel{\mathrm{˙}}{\mathit{\lambda }}}_{\mathrm{7}}& ={\mathit{\lambda }}_{\mathrm{8}}\\ {\stackrel{\mathrm{˙}}{\mathit{\lambda }}}_{\mathrm{8}}& =-\mathrm{2}{A}_{\mathrm{f}z}{\mathit{\lambda }}_{\mathrm{8}}-{B}_{\mathrm{f}z}{\mathit{\lambda }}_{\mathrm{7}}-\left(\mathrm{sin}{\mathit{\alpha }}_{\mathrm{n}}\\ & -\mathit{\mu }{\overline{f}}_{\mathit{\mu }}\mathrm{cos}{\mathit{\alpha }}_{\mathrm{n}}\right){B}_{\mathrm{df}}{f}_{\mathrm{b}}\left({\mathit{\lambda }}_{\mathrm{9}}\right)\\ & -\mathrm{2}\left(\mathrm{sin}{\mathit{\alpha }}_{\mathrm{n}}-\mathit{\mu }{\overline{f}}_{\mathit{\mu }}\mathrm{cos}{\mathit{\alpha }}_{\mathrm{n}}\right){A}_{\mathrm{df}}{\mathit{\lambda }}_{\mathrm{10}}\\ {\stackrel{\mathrm{˙}}{\mathit{\lambda }}}_{\mathrm{9}}& ={\mathit{\lambda }}_{\mathrm{10}}\\ {\stackrel{\mathrm{˙}}{\mathit{\lambda }}}_{\mathrm{10}}& =\mathrm{cos}{\mathit{\alpha }}_{\mathrm{n}}{\stackrel{\mathrm{˙}}{\mathit{\lambda }}}_{\mathrm{2}}+\mathrm{sin}{\mathit{\alpha }}_{\mathrm{n}}{\stackrel{\mathrm{˙}}{\mathit{\lambda }}}_{\mathrm{4}}\\ & -\mathrm{cos}{\mathit{\alpha }}_{\mathrm{n}}{\stackrel{\mathrm{˙}}{\mathit{\lambda }}}_{\mathrm{6}}-\mathrm{sin}{\mathit{\alpha }}_{\mathrm{n}}{\stackrel{\mathrm{˙}}{\mathit{\lambda }}}_{\mathrm{8}}-{B}_{\mathrm{d}e}\mathrm{cos}{\mathit{\alpha }}_{\mathrm{n}}{f}_{\mathrm{b}}\left({\mathit{\lambda }}_{\mathrm{9}}\right)\\ & -\mathrm{2}{A}_{\mathrm{d}e}\mathrm{cos}{\mathit{\alpha }}_{\mathrm{n}}{\mathit{\lambda }}_{\mathrm{10}}-\left(C+D\right)\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 Zc=17, ${Z}_{\mathrm{f}}=\mathrm{56}/\mathrm{68}/\mathrm{80}$. Modulus m=2mm. Pressure Angle αn=20. Input torque Tc=100N m. The dimensionless parameters of the vibration differential equation are assumed as: ${A}_{\mathrm{c}x}={A}_{\mathrm{c}z}=\mathrm{0.0261}$, Adc=0.093, ${B}_{\mathrm{c}x}={B}_{\mathrm{c}z}=\mathrm{0.0809}$, ${B}_{\mathrm{f}x}={B}_{\mathrm{f}z}=\mathrm{0.0071}$, Bdc=0.5509, ${A}_{\mathrm{f}x}={A}_{\mathrm{f}z}=\mathrm{3.2383}×{\mathrm{10}}^{-\mathrm{4}}$, ${A}_{\mathrm{df}}=\mathrm{3.8986}×{\mathrm{10}}^{-\mathrm{4}}$, Bdf=0.0231, Bde=0.8806, Ade=0.0149, $C=\mathrm{3.3982}×{\mathrm{10}}^{-\mathrm{5}}$, and $D=\mathrm{5.1342}×{\mathrm{10}}^{\mathrm{6}}×\mathrm{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 3Displacement, 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 4Displacement, phase, and Poincare diagrams of 1st stage face gear pair.

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

Figure 6Displacement, 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 7Displacement, phase, and Poincare diagrams of friction coefficient 0.1.

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

Figure 9Displacement, 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−1 respectively. The integrated meshing stiffness of the corresponding face gears pairs were 9.852×107N m−1, 5.7617×107N m−1, and 4.5731×107N m−1. The displacement curves, phase diagrams, and Poincare diagrams of the linear dynamics are shown in Figs. 10, 11 and 12.

Figure 10Displacement, phase, and Poincare diagrams of comprehensive meshing stiffness 9.852×7N m−1.

Figure 11Displacement, phase, and Poincare diagrams of comprehensive meshing stiffness 5.7617×107N m−1.

Figure 12Displacement, phase, and Poincare diagrams of comprehensive meshing stiffness 4.5731×107N m−1.

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−1, other parameters remained unchanged, the loading speed was set to 300, 500, 1000, and 3000 r min−1, 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−1.

## 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−1, the loading speed was 300 r min−1, 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 13Displacement, 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}$, Adc=0.093, ${B}_{\mathrm{c}x}={B}_{\mathrm{c}z}=\mathrm{0.0809}$, ${B}_{\mathrm{f}x}={B}_{\mathrm{f}z}=\mathrm{0.0071}$, Bdc=0.5509, ${A}_{\mathrm{f}x}={A}_{\mathrm{f}z}=\mathrm{3.2383}×{\mathrm{10}}^{-\mathrm{4}}$, ${A}_{\mathrm{df}}=\mathrm{3.8986}×{\mathrm{10}}^{-\mathrm{4}}$, Bdf=0.0231, Bde=0.8806, Ade=0.0149, $C=\mathrm{3.3982}×{\mathrm{10}}^{-\mathrm{5}}$, and $D=\mathrm{5.1342}×{\mathrm{10}}^{\mathrm{6}}×\mathrm{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−1, the nonlinear characteristics of face gears were analyzed, as shown in Fig. 14.

Figure 14Displacement, 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 15Displacement, phase, and Poincare diagrams of 1st stage face gear pair nonlinear system.

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

Figure 17Displacement, 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 18Displacement, 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−1, the loading speed was 300 r min−1, 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 19Displacement, 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−1, 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 20Displacement, phase, and Poincare diagrams of nonlinear system with excitation frequency 0.1648.

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

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

Figure 23Displacement, 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.

Data availability
Data availability.

Author contributions
Author contributions.

XC developed the overall concept of the paper, conducted the numerical analyses and wrote the majority of the paper. QH supervised the process and contributed in structuring the research, reviewing the results and reviewing the paper. ZX assisted in the experimental design and verified the suggested model. CZ checked the writing language.

Competing interests
Competing interests.

The authors declare that they have no conflict of interest.

Acknowledgements
Acknowledgements.

The authors would like to thank Xiaofeng Liang and Yao Ding for free help in English professional writing. This work was supported by the National Natural Science Foundation of China under Grant No. 51575191, and Guangdong Provincial Enterprises Key Laboratory Foundation about Mid-High-end Industrial Robot Technology under Grant No. 2018B030323027.

Financial support
Financial support.

This research has been supported by the National Natural Science Foundation of China (grant no. 51575191) and the Guangdong Provincial Enterprises Key Laboratory Foundation about Mid-High-end Industrial Robot Technology (grant no. 2018B030323027).

Review statement
Review statement.

This paper was edited by Guimin Chen and reviewed by Luca Bruzzone and one anonymous referee.

References

Atanasovska, I.: Multi-body contact in non-linear dynamics of real mechanical systems, X International Conference on Structural Dynamics (Eurodyn 2017), 199, 510–515, 2017.

Atanasovska, I. and Hedrih, K.: A New Collision Model for Analysing the Vibro-Impact of Spur Gears, T. Famena, 42, 1–13, 2018.

Bouslema, M., Frikha, A., Abdennadhar, M., Fakhfakh, T., Nasri, R., and Haddar, M.: Effects of modal truncation and condensation methods on the Frequency Response Function of a stage reducer connected by rigid coupling to a planetary gear system, Cr. Mecanique, 345, 807–823, 2017.

Cai, Z. Q. and Lin, C.: Dynamic Model and Analysis of Nonlinear Vibration Characteristic of a Curve-Face Gear Drive, Stroj. Vestn.-J. Mech. E, 63, 161–170, 2017.

Chowdhury, S. and Yedavalli, R. K.: Dynamics of low speed geared shaft systems mounted on rigid bearings, Mech. Mach. Theory, 112, 123–144, 2017.

Dadon, I., Koren, N., Klein, R., and Bortman, J.: A Step Toward Fault Type and Severity Characterization in Spur Gears, J, Mech, Design, 141, 8–30, 2019.

Dong, H. and Hu, Y. H.: Dynamic Load-Sharing Characteristic Analysis of Face Gear Power-Split Gear System Based on Tooth Contact Characteristics, Aip Conf. Proc., 1955, 030028-1–030028-5, 2018.

Fernandez-Del-Rincon, A., Garcia, P., Diez-Ibarbia, A., De-Juan, A., Iglesias, M., and Viadero, F.: Enhanced model of gear transmission dynamics for condition monitoring applications: Effects of torque, friction and bearing clearance, Mech. Syst. Signal Pr., 85, 445–467, 2017.

Hmida, A., Hammami, A., Chaari, F., Khabou, M. T., and Haddar, M.: Modal Analysis of Spur Gearbox with an Elastic Coupling, Advances in Acoustics and Vibration, 5, 153–163, 2017.

Hu, Z. H., Tang, J. Y., Chen, S. Y., and Sheng, Z. H.: Coupled translation-rotation vibration and dynamic analysis of face geared rotor system, J. Sound Vib., 351, 282–298, 2015.

Lin, C., Liu, Y., and Gu, S. J.: Analysis of nonlinear twisting vibration characteristics of orthogonal curve-face gear drive, J. Braz. Soc. Mech. Sci., 37, 1499–1505, 2015.

Liu, B., Zhao, J. G., Qian, J. H., and Wang, X. Y.: Nonlinear Vibration Analysis of Helical Face-Gear Transmission System with Multi-Factor, 2016 3rd International Conference on Mechanics and Mechatronics Research (Icmmr 2016), 77, 07010-16, 2016.

Liu, D. W., Gu, D. D., and Liu, Z. J.: Coupled vibration modeling and dynamic characteristics of noncircular face gear drive system with time-varying instantaneous center excitation, P. I. Mech. Eng. C-J. Mec., 233, 4947–4959, 2019.

Peng, M. and DeSmidt, H. A.: Torsional Stability of a Face-Gear Drive System, J. Am. Helicopter Soc., 60, 1–11, 2015.

Ren, F., Luo, G. F., Shi, G. Q., Wu, X. L., and Wang, N.: Influence of manufacturing errors on dynamic floating characteristics for herringbone planetary gears, Nonlinear Dynam., 93, 361–372, 2018.

Sakaridis, E., Spitas, V., and Spitas, C.: Non-linear modeling of gear drive dynamics incorporating intermittent tooth contact analysis and tooth eigenvibrations, Mech. Mach. Theory, 136, 307–333, 2019.

Saxena, A., Chouksey, M., and Parey, A.: Measurement of FRFs of coupled geared rotor system and the development of an accurate finite element model, Mech. Mach. Theory., 123, 66–75, 2018.

Tatar, A., Schwingshackl, C. W., and Friswell, M. I.: Dynamic behaviour of three-dimensional planetary geared rotor systems, Mech. Mach. Theory, 134, 39–56, 2019.

Wu, C. B., Cheng, J. H., and Zhang, G.: Study on Nonlinear Vibration of Planetary Gear System with NN – type Small Tooth Difference, 2017 32nd Youth Academic Annual Conference of Chinese Association of Automation (Yac), 2017, 645–649, 2017.

Xiao, Z. M., Wu, X., and Chen, Q.: Dynamic Modeling and Analysis of Multi-Stage Planetary Gears Coupled with Central Bearings and Housing, T. Can. Soc. Mech. Eng., 40, 1007–1018, 2016.

Yu, W. and Mechefske, C. K.: A New Model for the Single Mesh Stiffness Calculation of Helical Gears Using the Slicing Principle, Ijst-T. Mech. Eng., 43, 503–515, 2019.

Zajicek, M. and Dupal, J.: Analytical solution of spur gear mesh using linear model, Mech. Mach. Theory, 118, 154–167, 2017.

Zhang, A. Q., Wei, J., Qin, D. T., and Hou, S. S.: Analytical coupling characterization of multi-stage planetary gear free vibration considering flexible structure, J. Vibroeng., 19, 3994–4008, 2017.

Zhang, L., Wang, Y., Wu, K., Sheng, R. Y., and Huang, Q. L.: Dynamic modeling and vibration characteristics of a two-stage closed-form planetary gear train, Mech. Mach. Theory, 97, 12–28, 2016.

Zhang, L. N., Wang, Y., Wu, K., Sheng, R. Y., and Huang, Q. L.: Modal properties of a multi-stage power-split planetary gear set, Adv. Mech. Eng., 9, 1–10, 2017, 2017.