Establishment and analysis of nonlinear frequency response model of planetary gear transmission system

Based on the lumped parameter theory, a nonlinear bending torsion coupling dynamic model of planetary gear transmission system was established by considering the backlash, support clearance, time-varying meshing stiffness, meshing damping, transmission error and external periodic excitation. The model was solved by the Runge–Kutta method, the dynamic response was analyzed by a time domain diagram and phase diagram, and the nonlinear vibration characteristics were studied by the response curve of the speed vibration displacement. The vibration test of the planetary gearbox was carried out to verify the correctness of frequency domain response characteristics. The results show that the vibration response in the planetary gear system changes from a multiple periodic response to a single periodic response with the increase in input speed. Under the action of the backlash, time-varying meshing stiffness and meshing damping, the speed vibration displacement response curves of internal and external meshing pairs appear to form a nonlinear jump phenomenon and have a unilateral impact area, and the system presents nonlinear characteristics. The nonlinear vibration of the system can be effectively suppressed by decreasing the mesh stiffness or increasing the mesh resistance, while the vibration response displacement of the system increases by increasing the external exciting force, and the nonlinear characteristics of the system remain basically unchanged. The backlash is the main factor affecting the nonlinear frequency response of the system, but it can restrain the resonance of the system in a certain range. The spectrum characteristics of the vibration displacement signal of the planetary gearbox at different speeds are similar to the simulation results, which proves the validity of the simulation analysis model and the simulation results. It can provide a theoretical basis for the system vibration and noise reduction and a dynamic structural stability design optimization.

1 Introduction

Planetary gear transmission has been widely used in aircraft and vehicles because of its high load and large transmission ratio. However, due to its prominent nonlinear vibration problem, it is very prone to failure phenomena such as tooth fatigue pitting, tooth root crack and tooth or shaft fracture, which affects the operation accuracy, transmission efficiency and service life of the equipment.

Many scholars have done a lot of research on the nonlinear vibration of the planetary transmission system. Sun and Hu (2003) solved the bending and torsion model of planetary gear, considering the transverse support of solar shaft by the harmonic balance method, and verified the analysis results with a numerical simulation using the Runge–Kutta method. Tang et al. (2020) used the numerical integration method to analyze the effects of input speed, backlash, support clearance and damping ratio on the nonlinear dynamic characteristics of a planetary transmission system. Bahk and Parker (2011) solved and analyzed the vibration response and variation law of the system through perturbation analysis within the meaningful meshing frequency range and found that the external torque cannot inhibit the occurrence of the tooth falling off. Li et al. (2014) solved the multi-stage planetary dynamics model, based on the variable step size Gill numerical integration method, and analyzed the bifurcation characteristics of the system using the dynamic graphics tool. Li et al. (2016) considered the dynamic and static clearance and the frequency response characteristics of the gear system, and the influence of dynamic internal and external excitation parameters on the dynamic characteristics of the gear system were analyzed by using the analytical method. Li et al. (2015), in the study of the gear transmission system in a shearer, expressed and estimated the tooth clearance in the system by fractal expression, and the vibration characteristics of the system are studied theoretically and experimentally. Li et al. (2020) analyzed the theoretical fault characteristic frequency, combined the fault characteristic frequency of each component and diagnosed each combined fault based on Fourier spectrum analysis. Xiao et al. (2020) used the improved energy method to calculate the time-varying meshing stiffness of normal and cracked gears. Zhou et al. (2016) established a translational torsional coupling dynamic model, considering the time-varying meshing stiffness, comprehensive meshing error and backlash. Xiang et al. (2020a) analyzed the vibration response and bifurcation characteristics of the two-stage planetary gear train transmission system of a wind turbine gearbox. Jian et al. (2021) comprehensively considered the influence of the modification coefficient and time-varying meshing stiffness; based on the dynamic theory, the dynamic model of planetary gear system was established, and the thermal electrohydrodynamic lubrication characteristics of the modified gear system under vibration were analyzed. Chen et al. (2021) established the dynamic model of the planetary gearbox, considering the clearance of the planetary gear, sun gear and carrier bearing. Liu et al. (2020) studied and analyzed the influence of the flexibility of the ring gear and the support stiffness of the ring gear and the sun gear on the vibration characteristics of the planetary gear train. Wang et al. (2020) proposed a method to analyze the vibration and stress characteristics of the ring gear and analyzed the vibration and stress characteristics of the thin-walled ring gear. Li et al. (2020) established a multi-body dynamic model of the ring gear elastic support planetary gear system and discussed the influence of ring gear radial support stiffness on system vibration. Wang et al. (2021) proposed a vibration data analysis method based on the time synchronization average of enhanced planetary gears. Zhang et al. (2021) carried out a free modal test on the planetary transmission mechanism with the free modal test method and carried out a modal analysis with the finite element method. Xiang et al. (2020b) established a nonlinear model of a multi-stage gear transmission system, considering meshing stiffness, gear comprehensive error and backlash, and analyzed the nonlinear dynamic response of the system with excitation frequency and support stiffness as bifurcation parameters. Saghafi and Farshidianfar (2016) used the Melnikov method to give an analytical method to eliminate the chaos of the gear system under external control excitation. Yang et al. (2019) established a 3 degrees of freedom spur gear pair model, considering the nonlinearity of the backlash and bearing clearance, and studied the influence of tooth crack on the vibration characteristics of nonlinear spur gear system. Liu et al. (2017) studied the dynamic response of spur gear pair system and the interaction between the bearing clearance and backlash. Wang et al. (2018) obtained the bifurcation diagram and a Lyapunov exponent curve of the spur gear transmission system with the numerical method and analyzed the mechanism of the chaotic evolution of the gear transmission system. Fan et al. (2020) established a rigid flexible coupling dynamic model of the planetary gear transmission system, considering the flexibility of the inner ring gear and sun shaft, and calculated the dynamic response of the system with a numerical algorithm. Wang (2019) simulated the wear failure of the gear transmission system through the change of gear backlash and studied the bifurcation diagram of the system when the backlash increased.

To sum up, at present, the research on the nonlinear characteristics of the planetary transmission system mainly focuses on the bifurcation characteristics and fault analysis, while the research on the nonlinear frequency response characteristics of planetary gear system is less. In order to study the frequency response characteristics of the nonlinear system of the planetary gearbox under the action of multiple clearances, based on the lumped parameter theory and comprehensively considering the backlash, a nonlinear dynamic model of the bending torsion coupling of the planetary gear transmission system is established, including time-varying meshing stiffness, meshing synthesis error, external periodic excitation and other nonlinear factors. The model is solved by a numerical method, and the nonlinear dynamic characteristics and frequency response characteristics of the system are analyzed by the time domain diagram and phase diagram.

2 Establishment of nonlinear frequency response model

Using the lumped parameter method, the bending torsion coupling nonlinear vibration model of planetary gear transmission system, as shown in Fig. 1, is established.

When establishing the dynamic model, three coordinate systems, namely OXYZ, Oxyz and O_nx_ny_nz_n, are defined based on the sun gear, planet carrier and central axis of the planet gear. OXYZ is the fixed coordinate system, and Oxyz and O_nx_ny_nz_n are the rotational speed. ω_c is a moving coordinate system rotating around the axis O, and their coordinate axes are parallel to each other. The tooth surface friction is not considered in the bending torsion coupling nonlinear vibration model of the system. The inner ring gear is kept fixed. The influence of the input shaft of the sun gear and the output shaft of the planet carrier will temporarily not be considered. k_bjx, c_bjx and b_bjx represent support stiffness, support damping and the half-support clearance of the x axis. k_bjy, c_bjy and b_bjy represent the support stiffness, support damping and half-support clearance of the y axis. j = s, p and c represent the sun gear, planet gear and planet carrier, respectively. k_spi, c_spi, b_spi and e_spi represent the meshing stiffness, meshing damping, backlash and comprehensive transmission error of the meshing direction of the planet gear and the sun gear. k_rpi, c_rpi, b_rpi and e_rpi indicate the meshing stiffness, meshing damping, backlash and comprehensive transmission error of the meshing direction of planet gear and inner ring gear. Here, i=1, 2 and 3 denote the ith meshing pair.

Figure 1 Bending torsion coupling the nonlinear vibration model of the planetary gear transmission system.

3 Establishment of nonlinear frequency response differential equation

The model has 3N+6 degrees of freedom, where N is the number of the planet gear. The generalized displacement array X can be expressed as follows:

$$\begin{array}{}\text{(1)}& X={\left[x_{\mathrm{s}},y_{\mathrm{s}},{\mathit{\theta }}_{\mathrm{s}},x_{\mathrm{p}i},y_{\mathrm{p}i},{\mathit{\theta }}_{\mathrm{p}i},x_{\mathrm{c}},y_{\mathrm{c}},{\mathit{\theta }}_{\mathrm{c}}\right]}^{\mathrm{T}}.\end{array}$$

Here, x_j, y_j and θ_j are the vibration displacement, where j is s, pi and c.

The differential equations of the sun gear, the planet gear, and the planet carrier can be expressed as follows:

$$\begin{array}{}\text{(2)}& {\displaystyle }\left\{\begin{array}{l}{m}_{\mathrm{s}}{\ddot{x}}_{\mathrm{s}}=\sum _{i=\mathrm{1}}^{\mathrm{3}}\left[k_{\mathrm{sp}i}f\right({\mathit{\delta }}_{\mathrm{sp}i})+c_{\mathrm{sp}i}{\dot{\mathit{\delta }}}_{\mathrm{sp}i}]\sin \left({\mathit{\phi }}_{\mathrm{sp}i}\right)\\ -\left[k_{\mathrm{bs}x}f\right(x_{\mathrm{s}})+c_{\mathrm{bs}x}{\dot{x}}_{\mathrm{s}}]\\ m_{\mathrm{s}}{\ddot{y}}_{\mathrm{s}}=\sum _{i=\mathrm{1}}^{\mathrm{3}}\left[k_{\mathrm{sp}i}f\right({\mathit{\delta }}_{\mathrm{sp}i})+c_{\mathrm{sp}i}{\dot{\mathit{\delta }}}_{\mathrm{sp}i}]\cos \left({\mathit{\phi }}_{\mathrm{sp}i}\right)\\ -\left[k_{\mathrm{bs}y}f\right(y_{\mathrm{s}})+c_{\mathrm{bs}y}{\dot{y}}_{\mathrm{s}}]\\ I_{\mathrm{s}}{\ddot{\mathit{\theta }}}_{\mathrm{s}}=-\sum _{i=\mathrm{1}}^{\mathrm{3}}\left[k_{\mathrm{sp}i}f\right({\mathit{\delta }}_{\mathrm{sp}i})+c_{\mathrm{sp}i}{\dot{\mathit{\delta }}}_{\mathrm{sp}i}]r_{\mathrm{bs}}+T_{\mathrm{1}}\end{array}\right.\text{(3)}& {\displaystyle }\left\{\begin{array}{l}{m}_{\mathrm{p}}{\ddot{x}}_{\mathrm{p}i}=-\left[\right(k_{\mathrm{sp}i}f\left({\mathit{\delta }}_{\mathrm{sp}i}\right)+c_{\mathrm{sp}i}{\dot{\mathit{\delta }}}_{\mathrm{sp}i}\left)\sin \right({\mathit{\phi }}_{\mathrm{sp}i})\\ +\left(k_{\mathrm{rp}i}f\right({\mathit{\delta }}_{\mathrm{rp}i})+c_{\mathrm{rp}i}{\dot{\mathit{\delta }}}_{\mathrm{rp}i})\sin \left({\mathit{\phi }}_{\mathrm{rp}i}\right)]\\ -\left[k_{\mathrm{bp}x}f\right({\mathit{\delta }}_{\mathrm{cp}ix})+c_{\mathrm{bp}x}{\dot{\mathit{\delta }}}_{\mathrm{cp}ix}]\\ m_{\mathrm{p}}{\ddot{y}}_{\mathrm{p}i}=-\left[\right(k_{\mathrm{sp}i}f\left({\mathit{\delta }}_{\mathrm{sp}i}\right)+c_{\mathrm{sp}i}{\dot{\mathit{\delta }}}_{\mathrm{sp}i}\left)\cos \right({\mathit{\phi }}_{\mathrm{sp}i})\\ +\left(k_{\mathrm{rp}i}f\right({\mathit{\delta }}_{\mathrm{rp}i})+c_{\mathrm{rp}i}{\dot{\mathit{\delta }}}_{\mathrm{rp}i})\cos \left({\mathit{\phi }}_{\mathrm{rp}i}\right)]\\ -\left[k_{\mathrm{bp}y}f\right({\mathit{\delta }}_{\mathrm{cp}iy})+c_{\mathrm{bp}y}{\dot{\mathit{\delta }}}_{\mathrm{cp}iy}]\\ I_{\mathrm{p}}{\ddot{\mathit{\theta }}}_{\mathrm{p}i}=\left[k_{\mathrm{sp}i}f\right({\mathit{\delta }}_{\mathrm{sp}i})+c_{\mathrm{sp}i}{\dot{\mathit{\delta }}}_{\mathrm{sp}i}-k_{\mathrm{rp}i}f({\mathit{\delta }}_{\mathrm{rp}i})\\ -c_{\mathrm{rp}i}{\dot{\mathit{\delta }}}_{\mathrm{rp}i}]r_{\mathrm{bp}}\end{array}\right.\text{(4)}& {\displaystyle }\left\{\begin{array}{l}{m}_{\mathrm{c}}{\ddot{x}}_{\mathrm{c}}=\sum _{i=\mathrm{1}}^{\mathrm{3}}\left[k_{\mathrm{bp}x}f\right({\mathit{\delta }}_{\mathrm{cp}ix})+c_{\mathrm{bp}x}{\dot{\mathit{\delta }}}_{\mathrm{cp}ix}]\\ -[k_{\mathrm{bc}x}{x}_{\mathrm{c}}+c_{\mathrm{bc}x}{\dot{x}}_{\mathrm{c}}]\\ m_{\mathrm{c}}{\ddot{y}}_{\mathrm{c}}=\sum _{i=\mathrm{1}}^{\mathrm{3}}\left[k_{\mathrm{bp}y}f\right({\mathit{\delta }}_{\mathrm{cp}iy})+c_{\mathrm{bp}y}{\dot{\mathit{\delta }}}_{\mathrm{cp}iy}]\\ -[k_{\mathrm{bc}y}{x}_{\mathrm{c}}+c_{\mathrm{bc}y}{\dot{x}}_{\mathrm{c}}]\\ I_{\mathrm{c}}{\ddot{\mathit{\theta }}}_{\mathrm{c}}=-\sum _{i=\mathrm{1}}^{\mathrm{3}}\left[k_{\mathrm{bp}x}f\right({\mathit{\delta }}_{\mathrm{cp}ix})+c_{\mathrm{bp}x}{\dot{\mathit{\delta }}}_{\mathrm{cp}ix}]r_{\mathrm{bc}}\sin \left({\mathit{\phi }}_{\mathrm{p}i}\right)\\ +\left[k_{\mathrm{bp}y}f\right({\mathit{\delta }}_{\mathrm{cp}iy})+c_{\mathrm{bp}y}{\dot{\mathit{\delta }}}_{\mathrm{cp}iy}]r_{\mathrm{bc}}\cos \left({\mathit{\phi }}_{\mathrm{p}i}\right)+T_{\mathrm{2}}.\end{array}\right.\end{array}$$

Here, T₁ is the input torque. T₂ is the output torque. φ_pi is the position angle of the ith planet gear. m_j is the component mass. I_j is the moment of inertia of the component. r_bj is the radius of the base circle. j is s, p and c. φ_spi and φ_rpi are the angles between the y axis of the ith planet gear and its outer meshing plane and inner meshing plane, respectively. δ_spi and δ_rpi are the vibration displacement on the outer meshing line and the inner meshing line, respectively. δ_cpix and δ_cpiy are the vibration displacement of the ith planet gear relative to the planet carrier in each axis. i=1, 2, 3. f(δ) is the backlash nonlinear function.

The relative displacement in the direction of external meshing pair is defined as follows:

$$\begin{array}{}\text{(5)}& \begin{array}{rl}{\mathit{\delta }}_{\mathrm{sp}i}& =\left[x_{\mathrm{s}}\sin \right({\mathit{\phi }}_{\mathrm{sp}i})-x_{\mathrm{p}i}\sin ({\mathit{\phi }}_{\mathrm{p}i})+y_{\mathrm{s}}\cos ({\mathit{\phi }}_{\mathrm{sp}i})\\ & -y_{\mathrm{p}i}\cos \left({\mathit{\phi }}_{\mathrm{p}i}\right)+(r_{\mathrm{s}}{\mathit{\theta }}_{\mathrm{s}}+r_{\mathrm{p}i}{\mathit{\theta }}_{\mathrm{p}i})]-e_{\mathrm{sp}i}.\end{array}\end{array}$$

The relative displacement in the direction of internal meshing pair is defined as follows:

$$\begin{array}{}\text{(6)}& {\mathit{\delta }}_{\mathrm{rp}i}=\left(x_{\mathrm{p}i}\sin \right({\mathit{\phi }}_{\mathrm{p}i})-y_{\mathrm{p}i}\cos ({\mathit{\phi }}_{\mathrm{p}i})-r_{\mathrm{p}i}{\mathit{\theta }}_{\mathrm{p}i})-e_{\mathrm{rp}i}.\end{array}$$

The ith planet gear is on the x and y axis, relative to the planet carrier, as follows:

$$\begin{array}{}\text{(7)}& {\displaystyle }{\mathit{\delta }}_{\mathrm{cp}ix}=x_{\mathrm{p}i}-x_{\mathrm{c}}-r_{\mathrm{c}}{\mathit{\theta }}_{\mathrm{c}}\sin \left({\mathit{\phi }}_{\mathrm{p}i}\right)\text{(8)}& {\displaystyle }{\mathit{\delta }}_{\mathrm{cp}iy}=y_{\mathrm{p}i}-y_{\mathrm{c}}-r_{\mathrm{c}}{\mathit{\theta }}_{\mathrm{c}}\sin \left({\mathit{\phi }}_{\mathrm{p}i}\right).\end{array}$$

The system dynamics equation is dimensionless, and the nominal scale of $b_{\mathrm{c}}=\mathrm{10}\times {\mathrm{10}}^{-\mathrm{6}}$ and τ=ω_nt, of the order of $\mathrm{\Omega }=f_{\mathrm{m}}/{\mathit{\omega }}_{\mathrm{n}}$. ω_n is the excitation frequency. When the input speed n is given, the following expression can be obtained:

$$\begin{array}{}\text{(9)}& f_{\mathrm{m}}=f_{\mathrm{c}}{Z}_{\mathrm{r}}=(f_{\mathrm{s}}^{\left(\mathrm{r}\right)}-f_{\mathrm{c}})Z_{\mathrm{s}}.\end{array}$$

Here, Z_s and Z_r are the number of sun gear teeth and inner ring gear teeth, respectively, $f_{\mathrm{s}}^{\left(\mathrm{r}\right)}$ is the absolute rotation frequency of the sun gear, and f_c is the rotation frequency of planet carrier. The calculation formula is as follows:

$$\begin{array}{}\text{(10)}& {\displaystyle }{f}_{\mathrm{s}}^{\left(\mathrm{r}\right)}=n/\mathrm{60}\text{(11)}& {\displaystyle }{f}_{\mathrm{c}}=f_{\mathrm{s}}^{\left(\mathrm{r}\right)}/(\mathrm{1}-(-Z_{\mathrm{r}}/Z_{\mathrm{s}}\left)\right).\end{array}$$

The formula of ω_n is as follows:

$$\begin{array}{}\text{(12)}& {\mathit{\omega }}_{\mathrm{n}}=\sqrt{k_{\mathrm{sp}i}/(\mathrm{1}/m_{\mathrm{es}}+\mathrm{1}/m_{\mathrm{ep}}+\mathrm{1}/m_{\mathrm{s}}+\mathrm{1}/m_{\mathrm{p}})}.\end{array}$$

Here, m_es and m_ep are the equivalent masses of the sun gear and planet gear, respectively, $m_{\mathrm{es}}=I_{\mathrm{s}}/r_{\mathrm{bs}}$, and $m_{\mathrm{ep}}=I_{\mathrm{p}i}/r_{\mathrm{bp}}$.

Other dimensionless physical quantities are defined as follows:

$$\begin{array}{}\text{(13)}& \left\{\begin{array}{l}\stackrel{\mathrm{‾}}{X}\left(\mathit{\tau }\right)=X/b_{\mathrm{c}},\dot{\stackrel{\mathrm{‾}}{X}}\left(\mathit{\tau }\right)=\dot{X}/\left(b_{\mathrm{c}}{\mathit{\omega }}_{\mathrm{n}}\right),\ddot{\stackrel{\mathrm{‾}}{X}}\left(\mathit{\tau }\right)=\ddot{X}/\left(b_{\mathrm{c}}{\mathit{\omega }}_{\mathrm{n}}^{\mathrm{2}}\right),\\ \\ \stackrel{\mathrm{‾}}{e}\left(\mathit{\tau }\right)=e/b_{\mathrm{c}},\dot{\stackrel{\mathrm{‾}}{e}}\left(\mathit{\tau }\right)=e/\left(b_{\mathrm{c}}{\mathit{\omega }}_{\mathrm{n}}\right),\ddot{\stackrel{\mathrm{‾}}{e}}\left(\mathit{\tau }\right)=e/\left(b_{\mathrm{c}}{\mathit{\omega }}_{\mathrm{n}}^{\mathrm{2}}\right),\\ \\ \stackrel{\mathrm{‾}}{c}=c/{\mathit{\omega }}_{\mathrm{n}},\stackrel{\mathrm{‾}}{k}=k/{\mathit{\omega }}_{\mathrm{n}}^{\mathrm{2}},{\stackrel{\mathrm{‾}}{T}}_{\mathrm{1}}=T_{\mathrm{1}}/b_{\mathrm{c}}{\mathit{\omega }}_{\mathrm{n}}^{\mathrm{2}},{\stackrel{\mathrm{‾}}{T}}_{\mathrm{2}}=T_{\mathrm{2}}/b_{\mathrm{c}}{\mathit{\omega }}_{\mathrm{n}}^{\mathrm{2}},\\ \\ \stackrel{\mathrm{‾}}{b}=b/b_{\mathrm{c}},F_{\mathrm{1}}=\stackrel{\mathrm{‾}}{T}/m_{\mathrm{es}}{r}_{\mathrm{s}},F_{\mathrm{2}}={\stackrel{\mathrm{‾}}{T}}_{\mathrm{2}}/m_{\mathrm{ec}}{r}_{\mathrm{c}}.\end{array}\right.\end{array}$$

The external exciting force F₁ is expanded by Fourier series, and the order is the same as the stiffness, which can be expressed as follows:

$$\begin{array}{}\text{(14)}& F_{\mathrm{m}}\left(\stackrel{\mathrm{‾}}{t}\right)=F_{\mathrm{1}}+{\sum }_{l=\mathrm{1}}^L{f}_l\cos \left(l\mathrm{\Omega }\stackrel{\mathrm{‾}}{t}\right).\end{array}$$

Here, f_l is the lth harmonic coefficient of the input dynamic excitation ${\stackrel{\mathrm{‾}}{T}}_{\mathrm{1}}$ expanded by Fourier series.

By introducing the relative displacement and dimensionless physical quantity defined above into the differential equation, the dimensionless differential equation with $\stackrel{\mathrm{‾}}{X}\left(\mathit{\tau }\right)$ as the variable is obtained as follows:

$$\begin{array}{}\text{(15)}& \stackrel{\mathrm{‾}}{X}\left(\mathit{\tau }\right)={\left[{\stackrel{\mathrm{‾}}{x}}_{\mathrm{s}},{\stackrel{\mathrm{‾}}{y}}_{\mathrm{s}},{\stackrel{\mathrm{‾}}{x}}_{\mathrm{p}i},{\stackrel{\mathrm{‾}}{y}}_{\mathrm{p}i},{\stackrel{\mathrm{‾}}{\mathit{\delta }}}_{\mathrm{sp}i},{\stackrel{\mathrm{‾}}{\mathit{\delta }}}_{\mathrm{rp}i}\right]}^{\mathrm{T}}.\end{array}$$

The dimensionless differential equations of the sun gear and the planet gear are expressed as follows:

$$\begin{array}{}\text{(16)}& {\displaystyle }\left\{\begin{array}{l}{\ddot{\stackrel{\mathrm{‾}}{x}}}_{\mathrm{s}}=\sum _{i=\mathrm{1}}^{\mathrm{3}}[\frac{k_{\mathrm{sp}i}f\left({\stackrel{\mathrm{‾}}{\mathit{\delta }}}_{\mathrm{sp}i}\right)}{{\mathit{\omega }}_{\mathrm{n}}^{\mathrm{2}}{m}_{\mathrm{s}}}+\frac{c_{\mathrm{sp}i}{\dot{\stackrel{\mathrm{‾}}{\mathit{\delta }}}}_{\mathrm{sp}i}}{{\mathit{\omega }}_{\mathrm{n}}{m}_{\mathrm{s}}}]\sin \left({\mathit{\phi }}_{\mathrm{sp}i}\right)\\ -[\frac{k_{\mathrm{bs}x}f\left({\stackrel{\mathrm{‾}}{x}}_{\mathrm{s}}\right)}{{\mathit{\omega }}_{\mathrm{n}}^{\mathrm{2}}{m}_{\mathrm{s}}}+\frac{c_{\mathrm{bs}x}{\dot{\stackrel{\mathrm{‾}}{x}}}_{\mathrm{s}}}{{\mathit{\omega }}_{\mathrm{n}}{m}_{\mathrm{s}}}]\\ {\ddot{\stackrel{\mathrm{‾}}{y}}}_{\mathrm{s}}=\sum _{i=\mathrm{1}}^{\mathrm{3}}[\frac{k_{\mathrm{sp}i}f\left({\stackrel{\mathrm{‾}}{\mathit{\delta }}}_{\mathrm{sp}i}\right)}{{\mathit{\omega }}_{\mathrm{n}}^{\mathrm{2}}{m}_{\mathrm{s}}}+\frac{c_{\mathrm{sp}i}{\dot{\stackrel{\mathrm{‾}}{\mathit{\delta }}}}_{\mathrm{sp}i}}{{\mathit{\omega }}_{\mathrm{n}}{m}_{\mathrm{s}}}]\cos \left({\mathit{\phi }}_{\mathrm{sp}i}\right)\\ -[\frac{k_{\mathrm{bs}y}f\left({\stackrel{\mathrm{‾}}{y}}_{\mathrm{s}}\right)}{{\mathit{\omega }}_{\mathrm{n}}^{\mathrm{2}}{m}_{\mathrm{s}}}+\frac{c_{\mathrm{bs}y}{\dot{\stackrel{\mathrm{‾}}{y}}}_{\mathrm{s}}}{{\mathit{\omega }}_{\mathrm{n}}{m}_{\mathrm{s}}}]\end{array}\right.\text{(17)}& {\displaystyle }\left\{\begin{array}{l}{\ddot{\stackrel{\mathrm{‾}}{x}}}_{\mathrm{p}i}=-\left[\left(\frac{k_{\mathrm{sp}i}f\left({\stackrel{\mathrm{‾}}{\mathit{\delta }}}_{\mathrm{sp}i}\right)}{{\mathit{\omega }}_{\mathrm{n}}^{\mathrm{2}}{m}_{\mathrm{p}i}}+\frac{c_{\mathrm{sp}i}{\dot{\stackrel{\mathrm{‾}}{\mathit{\delta }}}}_{\mathrm{sp}i}}{{\mathit{\omega }}_{\mathrm{n}}{m}_{\mathrm{p}i}}\right)\sin \right({\mathit{\phi }}_{\mathrm{sp}i})\\ +\left(\frac{k_{\mathrm{rp}i}f\left({\stackrel{\mathrm{‾}}{\mathit{\delta }}}_{\mathrm{rp}i}\right)}{{\mathit{\omega }}_{\mathrm{n}}^{\mathrm{2}}{m}_{\mathrm{p}i}}+\frac{c_{\mathrm{rp}i}{\dot{\stackrel{\mathrm{‾}}{\mathit{\delta }}}}_{\mathrm{rp}i}}{{\mathit{\omega }}_{\mathrm{n}}{m}_{\mathrm{p}i}}\right)\sin \left({\mathit{\phi }}_{\mathrm{rp}i}\right)]\\ -[\frac{k_{\mathrm{bp}x}f\left({\stackrel{\mathrm{‾}}{\mathit{\delta }}}_{\mathrm{cp}ix}\right)}{{\mathit{\omega }}_{\mathrm{n}}^{\mathrm{2}}{m}_{\mathrm{p}i}}+\frac{c_{\mathrm{bp}x}{\dot{\stackrel{\mathrm{‾}}{\mathit{\delta }}}}_{\mathrm{cp}ix}}{{\mathit{\omega }}_{\mathrm{n}}{m}_{\mathrm{p}i}}]\\ {\ddot{\stackrel{\mathrm{‾}}{y}}}_{\mathrm{p}i}=-\left[\left(\frac{k_{\mathrm{sp}i}f\left({\stackrel{\mathrm{‾}}{\mathit{\delta }}}_{\mathrm{sp}i}\right)}{{\mathit{\omega }}_{\mathrm{n}}^{\mathrm{2}}{m}_{\mathrm{p}i}}+\frac{c_{\mathrm{sp}i}{\dot{\stackrel{\mathrm{‾}}{\mathit{\delta }}}}_{\mathrm{sp}i}}{{\mathit{\omega }}_{\mathrm{n}}{m}_{\mathrm{p}i}}\right)\cos \right({\mathit{\phi }}_{\mathrm{sp}i})\\ +\left(\frac{k_{\mathrm{rp}i}f\left({\stackrel{\mathrm{‾}}{\mathit{\delta }}}_{\mathrm{rp}i}\right)}{{\mathit{\omega }}_{\mathrm{n}}^{\mathrm{2}}{m}_{\mathrm{p}i}}+\frac{c_{\mathrm{rp}i}{\dot{\stackrel{\mathrm{‾}}{\mathit{\delta }}}}_{\mathrm{rp}i}}{{\mathit{\omega }}_{\mathrm{n}}{m}_{\mathrm{p}i}}\right)\cos \left({\mathit{\phi }}_{\mathrm{rp}i}\right)]\\ -[\frac{k_{\mathrm{bp}y}f\left({\stackrel{\mathrm{‾}}{\mathit{\delta }}}_{\mathrm{cp}iy}\right)}{{\mathit{\omega }}_{\mathrm{n}}^{\mathrm{2}}{m}_{\mathrm{p}i}}+\frac{c_{\mathrm{bp}y}{\dot{\stackrel{\mathrm{‾}}{\mathit{\delta }}}}_{\mathrm{cp}iy}}{{\mathit{\omega }}_{\mathrm{n}}{m}_{\mathrm{p}i}}].\end{array}\right.\end{array}$$

The dimensionless differential equation of the internal and external meshing pair can be expressed as follows:

$$\begin{array}{}\text{(18)}& {\displaystyle }\begin{array}{l}{\ddot{\stackrel{\mathrm{‾}}{\mathit{\delta }}}}_{\mathrm{rp}i}=\left(\right[(\frac{k_{\mathrm{sp}i}f\left({\stackrel{\mathrm{‾}}{\mathit{\delta }}}_{\mathrm{sp}i}\right)}{{\mathit{\omega }}_{\mathrm{n}}^{\mathrm{2}}{m}_{\mathrm{p}i}}+\frac{c_{\mathrm{sp}i}{\dot{\stackrel{\mathrm{‾}}{\mathit{\delta }}}}_{\mathrm{sp}i}}{{\mathit{\omega }}_{\mathrm{n}}{m}_{\mathrm{p}i}})\sin \left({\mathit{\phi }}_{\mathrm{sp}i}\right)\\ -(\frac{k_{\mathrm{rp}i}f\left({\stackrel{\mathrm{‾}}{\mathit{\delta }}}_{\mathrm{rp}i}\right)}{{\mathit{\omega }}_{\mathrm{n}}^{\mathrm{2}}{m}_{\mathrm{p}i}}+\frac{c_{\mathrm{rp}i}{\dot{\stackrel{\mathrm{‾}}{\mathit{\delta }}}}_{\mathrm{rp}i}}{{\mathit{\omega }}_{\mathrm{n}}{m}_{\mathrm{p}i}})\sin \left({\mathit{\phi }}_{\mathrm{rp}i}\right)]\\ +[\frac{k_{\mathrm{bp}x}f\left({\stackrel{\mathrm{‾}}{\mathit{\delta }}}_{\mathrm{cp}ix}\right)}{{\mathit{\omega }}_{\mathrm{n}}^{\mathrm{2}}{m}_{\mathrm{p}i}}+\frac{c_{\mathrm{bp}x}{\dot{\stackrel{\mathrm{‾}}{\mathit{\delta }}}}_{\mathrm{cp}ix}}{{\mathit{\omega }}_{\mathrm{n}}{m}_{\mathrm{p}i}}]\left)\sin \right({\mathit{\phi }}_{\mathrm{p}i})\\ -\left(\right[(\frac{k_{\mathrm{sp}i}f\left({\stackrel{\mathrm{‾}}{\mathit{\delta }}}_{\mathrm{sp}i}\right)}{{\mathit{\omega }}_{\mathrm{n}}^{\mathrm{2}}{m}_{\mathrm{p}i}}+\frac{c_{\mathrm{sp}i}{\dot{\stackrel{\mathrm{‾}}{\mathit{\delta }}}}_{\mathrm{sp}i}}{{\mathit{\omega }}_{\mathrm{n}}{m}_{\mathrm{p}i}})\cos \left({\mathit{\phi }}_{\mathrm{sp}i}\right)\\ -(\frac{k_{\mathrm{rp}i}f\left({\stackrel{\mathrm{‾}}{\mathit{\delta }}}_{\mathrm{rp}i}\right)}{{\mathit{\omega }}_{\mathrm{n}}^{\mathrm{2}}{m}_{\mathrm{p}i}}+\frac{c_{\mathrm{rp}i}{\dot{\stackrel{\mathrm{‾}}{\mathit{\delta }}}}_{\mathrm{rp}i}}{{\mathit{\omega }}_{\mathrm{n}}{m}_{\mathrm{p}i}})\cos \left({\mathit{\phi }}_{\mathrm{rp}i}\right)]\\ +[\frac{k_{\mathrm{bp}y}f\left({\stackrel{\mathrm{‾}}{\mathit{\delta }}}_{\mathrm{cp}iy}\right)}{{\mathit{\omega }}_{\mathrm{n}}^{\mathrm{2}}{m}_{\mathrm{p}i}}+\frac{c_{\mathrm{bp}y}{\dot{\stackrel{\mathrm{‾}}{\mathit{\delta }}}}_{\mathrm{cp}iy}}{{\mathit{\omega }}_{\mathrm{n}}{m}_{\mathrm{p}i}}]\left)\cos \right({\mathit{\phi }}_{\mathrm{p}i})\\ -[\frac{k_{\mathrm{sp}i}f\left({\stackrel{\mathrm{‾}}{\mathit{\delta }}}_{\mathrm{sp}i}\right)}{{\mathit{\omega }}_{\mathrm{n}}^{\mathrm{2}}{m}_{\mathrm{ep}}}+\frac{c_{\mathrm{sp}i}{\dot{\stackrel{\mathrm{‾}}{\mathit{\delta }}}}_{\mathrm{sp}i}}{{\mathit{\omega }}_{\mathrm{n}}{m}_{\mathrm{ep}}}-\frac{k_{\mathrm{rp}i}f\left({\stackrel{\mathrm{‾}}{\mathit{\delta }}}_{\mathrm{rp}i}\right)}{{\mathit{\omega }}_{\mathrm{n}}^{\mathrm{2}}{m}_{\mathrm{ep}}}-\frac{c_{\mathrm{rp}i}{\dot{\stackrel{\mathrm{‾}}{\mathit{\delta }}}}_{\mathrm{rp}i}}{{\mathit{\omega }}_{\mathrm{n}}{m}_{\mathrm{ep}}}]-{\ddot{\stackrel{\mathrm{‾}}{e}}}_{\mathrm{rp}i}\end{array}\text{(19)}& {\displaystyle }\begin{array}{l}{\ddot{\stackrel{\mathrm{‾}}{\mathit{\delta }}}}_{\mathrm{sp}i}=\left(\sum _{i=\mathrm{1}}^{\mathrm{3}}\right[\frac{k_{\mathrm{sp}i}f\left({\stackrel{\mathrm{‾}}{\mathit{\delta }}}_{\mathrm{sp}i}\right)}{{\mathit{\omega }}_{\mathrm{n}}^{\mathrm{2}}{m}_{\mathrm{s}}}+\frac{c_{\mathrm{sp}i}{\dot{\stackrel{\mathrm{‾}}{\mathit{\delta }}}}_{\mathrm{sp}i}}{{\mathit{\omega }}_{\mathrm{n}}{m}_{\mathrm{s}}}\left]\sin \right({\mathit{\phi }}_{\mathrm{sp}i})\\ -[\frac{k_{\mathrm{bs}x}f\left({\stackrel{\mathrm{‾}}{x}}_{\mathrm{s}}\right)}{{\mathit{\omega }}_{\mathrm{n}}^{\mathrm{2}}{m}_{\mathrm{s}}}+\frac{c_{\mathrm{bs}x}{\dot{\stackrel{\mathrm{‾}}{x}}}_{\mathrm{s}}}{{\mathit{\omega }}_{\mathrm{n}}{m}_{\mathrm{s}}}]\left)\sin \right({\mathit{\phi }}_{\mathrm{sp}i})\\ +\left(\right[(\frac{k_{\mathrm{sp}i}f\left({\stackrel{\mathrm{‾}}{\mathit{\delta }}}_{\mathrm{sp}i}\right)}{{\mathit{\omega }}_{\mathrm{n}}^{\mathrm{2}}{m}_{\mathrm{p}i}}+\frac{c_{\mathrm{sp}i}{\dot{\stackrel{\mathrm{‾}}{\mathit{\delta }}}}_{\mathrm{sp}i}}{{\mathit{\omega }}_{\mathrm{n}}{m}_{\mathrm{p}i}})\sin \left({\mathit{\phi }}_{\mathrm{sp}i}\right)\\ +(\frac{k_{\mathrm{rp}i}f\left({\stackrel{\mathrm{‾}}{\mathit{\delta }}}_{\mathrm{rp}i}\right)}{{\mathit{\omega }}_{\mathrm{n}}^{\mathrm{2}}{m}_{\mathrm{p}i}}+\frac{c_{\mathrm{rp}i}{\dot{\stackrel{\mathrm{‾}}{\mathit{\delta }}}}_{\mathrm{rp}i}}{{\mathit{\omega }}_{\mathrm{n}}{m}_{\mathrm{p}i}})\sin \left({\mathit{\phi }}_{\mathrm{rp}i}\right)]\\ +[\frac{k_{\mathrm{bp}x}f\left({\stackrel{\mathrm{‾}}{\mathit{\delta }}}_{\mathrm{cp}ix}\right)}{{\mathit{\omega }}_{\mathrm{n}}^{\mathrm{2}}{m}_{\mathrm{p}i}}+\frac{c_{\mathrm{bp}x}{\dot{\stackrel{\mathrm{‾}}{\mathit{\delta }}}}_{\mathrm{cp}ix}}{{\mathit{\omega }}_{\mathrm{n}}{m}_{\mathrm{p}i}}]\left)\sin \right({\mathit{\phi }}_{\mathrm{p}i})\\ +\left(\sum _{i=\mathrm{1}}^{\mathrm{3}}\right[\frac{k_{\mathrm{sp}i}f\left({\stackrel{\mathrm{‾}}{\mathit{\delta }}}_{\mathrm{sp}i}\right)}{{\mathit{\omega }}_{\mathrm{n}}^{\mathrm{2}}{m}_{\mathrm{s}}}+\frac{c_{\mathrm{sp}i}{\dot{\stackrel{\mathrm{‾}}{\mathit{\delta }}}}_{\mathrm{sp}i}}{{\mathit{\omega }}_{\mathrm{n}}{m}_{\mathrm{s}}}\left]\cos \right({\mathit{\phi }}_{\mathrm{sp}i})\\ -[\frac{k_{\mathrm{bs}y}f\left({\stackrel{\mathrm{‾}}{y}}_{\mathrm{s}}\right)}{{\mathit{\omega }}_{\mathrm{n}}^{\mathrm{2}}{m}_{\mathrm{s}}}+\frac{c_{\mathrm{bs}y}{\dot{\stackrel{\mathrm{‾}}{y}}}_{\mathrm{s}}}{{\mathit{\omega }}_{\mathrm{n}}{m}_{\mathrm{s}}}]\left)\cos \right({\mathit{\phi }}_{\mathrm{sp}i})\\ -\left(\right(\frac{k_{\mathrm{sp}i}f\left({\stackrel{\mathrm{‾}}{\mathit{\delta }}}_{\mathrm{sp}i}\right)}{{\mathit{\omega }}_{\mathrm{n}}^{\mathrm{2}}{m}_{\mathrm{p}i}}+\frac{c_{\mathrm{sp}i}{\dot{\stackrel{\mathrm{‾}}{\mathit{\delta }}}}_{\mathrm{sp}i}}{{\mathit{\omega }}_{\mathrm{n}}{m}_{\mathrm{p}i}}\left)\cos \right({\mathit{\phi }}_{\mathrm{sp}i})\\ +(\frac{k_{\mathrm{rp}i}f\left({\stackrel{\mathrm{‾}}{\mathit{\delta }}}_{\mathrm{rp}i}\right)}{{\mathit{\omega }}_{\mathrm{n}}^{\mathrm{2}}{m}_{\mathrm{p}i}}+\frac{c_{\mathrm{rp}i}{\dot{\stackrel{\mathrm{‾}}{\mathit{\delta }}}}_{\mathrm{rp}i}}{{\mathit{\omega }}_{\mathrm{n}}{m}_{\mathrm{p}i}})\cos \left({\mathit{\phi }}_{\mathrm{rp}i}\right)\\ +[\frac{k_{\mathrm{bp}y}f\left({\stackrel{\mathrm{‾}}{\mathit{\delta }}}_{\mathrm{cp}iy}\right)}{{\mathit{\omega }}_{\mathrm{n}}^{\mathrm{2}}{m}_{\mathrm{p}i}}+\frac{c_{\mathrm{bp}y}{\dot{\stackrel{\mathrm{‾}}{\mathit{\delta }}}}_{\mathrm{cp}iy}}{{\mathit{\omega }}_{\mathrm{n}}{m}_{\mathrm{p}i}}]\left)\cos \right({\mathit{\phi }}_{\mathrm{p}i})\\ -\left(\sum _{i=\mathrm{1}}^{\mathrm{3}}\right[\frac{k_{\mathrm{sp}i}f\left({\stackrel{\mathrm{‾}}{\mathit{\delta }}}_{\mathrm{sp}i}\right)}{{\mathit{\omega }}_{\mathrm{n}}^{\mathrm{2}}}+\frac{c_{\mathrm{sp}i}{\dot{\stackrel{\mathrm{‾}}{\mathit{\delta }}}}_{\mathrm{sp}i}}{{\mathit{\omega }}_{\mathrm{n}}}]-\frac{T_{\mathrm{1}}}{\mathrm{2}{r}_{\mathrm{bs}}{m}_{\mathrm{es}}{\mathit{\omega }}_{\mathrm{n}}^{\mathrm{2}}{b}_{\mathrm{c}}})\\ +[\frac{k_{\mathrm{sp}i}f\left({\stackrel{\mathrm{‾}}{\mathit{\delta }}}_{\mathrm{sp}i}\right)}{{\mathit{\omega }}_{\mathrm{n}}^{\mathrm{2}}{m}_{\mathrm{ep}}}+\frac{c_{\mathrm{sp}i}{\dot{\stackrel{\mathrm{‾}}{\mathit{\delta }}}}_{\mathrm{sp}i}}{{\mathit{\omega }}_{\mathrm{n}}{m}_{\mathrm{ep}}}-\frac{k_{\mathrm{rp}i}f\left({\stackrel{\mathrm{‾}}{\mathit{\delta }}}_{\mathrm{rp}i}\right)}{{\mathit{\omega }}_{\mathrm{n}}^{\mathrm{2}}{m}_{\mathrm{ep}}}-\frac{c_{\mathrm{rp}i}{\dot{\stackrel{\mathrm{‾}}{\mathit{\delta }}}}_{\mathrm{rp}i}}{{\mathit{\omega }}_{\mathrm{n}}{m}_{\mathrm{ep}}}]-{\ddot{\stackrel{\mathrm{‾}}{e}}}_{\mathrm{sp}i}.\end{array}\end{array}$$

4 Example of system frequency response characteristics and parameter influence analysis

4.1 Frequency response characteristics analysis

The basic parameters are shown in Table 1.

Table 1 Basic parameters of planetary system.

The time-varying meshing stiffness k_rpi and k_spi can be expressed with the following equations:

$$\begin{array}{}\text{(20)}& \begin{array}{rl}& k_{\mathrm{rp}i}\left(t\right)=k_{\mathrm{rm}}\left(\mathrm{1}+\sum _{b=\mathrm{1}}^L{k}_b\cos (b{f}_{\mathrm{m}}t+{\mathit{\lambda }}_b)\right)\\ & k_{\mathrm{sp}i}\left(t\right)=k_{\mathrm{sm}}\left(\mathrm{1}+\sum _{a=\mathrm{1}}^L{k}_a\cos (a{f}_{\mathrm{m}}t+{\mathit{\lambda }}_a)\right).\end{array}\end{array}$$

Here, k_sm and k_rm are the average meshing stiffness. k_a and k_b are the harmonic coefficients of the meshing stiffness. λ_a and λ_b are the phase angles of the meshing stiffness.

The comprehensive transmission errors e_rpi and e_spi can be expressed with the following equations:

$$\begin{array}{}\text{(21)}& \begin{array}{rl}& e_{\mathrm{rp}i}\left(t\right)=e_{\mathrm{r}}\left(\mathrm{1}+\sum _{b=\mathrm{1}}^L{e}_b\sin (b{f}_{\mathrm{m}}t+{\mathit{\phi }}_b)\right)\\ & e_{\mathrm{sp}i}\left(t\right)=e_{\mathrm{s}}\left(\mathrm{1}+\sum _{a=\mathrm{1}}^L{e}_a\sin (a{f}_{\mathrm{m}}t+{\mathit{\phi }}_a)\right).\end{array}\end{array}$$

Here, e_s and e_r are the fluctuating mean values of the meshing comprehensive errors. e_a and e_b are harmonic coefficients of the meshing comprehensive errors. φ_a and φ_b are the phase angles of the meshing comprehensive errors.

The nonlinear function of the meshing backlash and support clearance can be defined as follows:

$$\begin{array}{}\text{(22)}& \begin{array}{rl}& f\left({\mathit{\delta }}_{\mathrm{sp}i}\right)=\left\{\begin{array}{ll}{\mathit{\delta }}_{\mathrm{sp}i}-b_{\mathrm{sp}i}& {\mathit{\delta }}_{\mathrm{sp}i}>b_{\mathrm{sp}i}\\ \mathrm{0}& \left|{\mathit{\delta }}_{\mathrm{sp}i}\right|\le b_{\mathrm{sp}i}\\ {\mathit{\delta }}_{\mathrm{sp}i}+b_{\mathrm{sp}i}& {\mathit{\delta }}_{\mathrm{sp}i}<-b_{\mathrm{sp}i};\end{array}\right.\\ & f\left({\mathit{\delta }}_{\mathrm{rp}i}\right)=\left\{\begin{array}{ll}{\mathit{\delta }}_{\mathrm{rp}i}-b_{\mathrm{rp}i}& {\mathit{\delta }}_{\mathrm{rp}i}>b_{\mathrm{rp}i}\\ \mathrm{0}& \left|{\mathit{\delta }}_{\mathrm{rp}i}\right|\le b_{\mathrm{rp}i}\\ {\mathit{\delta }}_{\mathrm{rp}i}+b_{\mathrm{rp}i}& {\mathit{\delta }}_{\mathrm{rp}i}<-b_{\mathrm{rp}i}.\end{array}\right.\end{array}\end{array}$$

The damping coefficients can be expressed as follows:

$$\begin{array}{}\text{(23)}& \begin{array}{rl}& c_{\mathrm{sp}i}=\mathrm{2}{\mathit{\xi }}_{\mathrm{sp}i}\sqrt{{\displaystyle \frac{k_{\mathrm{sp}i}}{\mathrm{1}/m_{\mathrm{es}}+\mathrm{1}/m_{\mathrm{ep}}}}};\\ & c_{\mathrm{rp}i}=\mathrm{2}{\mathit{\xi }}_{\mathrm{rp}i}\sqrt{{\displaystyle \frac{k_{\mathrm{rp}i}}{\mathrm{1}/m_{\mathrm{er}}+\mathrm{1}/m_{\mathrm{ep}}}}}.\end{array}\end{array}$$

Here, ξ_rpi and ξ_spi are the meshing damping ratios of the ith internal and external meshing pairs.

The harmonic parameters of each excitation term in the external meshing pair can be expressed as follows: L=3, N=3, F₁=1, f₁=0.05, f₂=0.02, f₃=0.01, e₁=0.05, e₂=0.02, e₃=0.01, k₁=0.1, k₂=0.04, and k₃=0.05. The harmonic parameters of each excitation term in the internal meshing pair can be expressed as follows: L=3, N=3, F₁=1, f₁=0.05, f₂=0.02, f₃=0.01, e₁=0.05, e₂=0.02, e₃=0.01, k₁=0.05, k₂=0.02, and k₃=0.01.

Table 2 Dynamic parameters of planetary system.

The average meshing stiffness of the internal and external gear pairs is calculated according to the China national product standard (GB/T 3480-2019). The calculation results are shown in Table 2 and are verified by the Ishikawa method (Mao et al., 2021). The calculated meshing stiffness is expanded by a Fourier series to obtain the internal and external time-varying meshing stiffness curve, as shown in Fig. 2.

Figure 2 Time-varying meshing stiffness of the planetary system.

The Runge–Kutta method is used to solve the differential equation to obtain the vibration response of the internal and external meshing pairs of the planetary system at different speeds, as shown in Figs. 3 and 4.

Figure 3 Vibration response of the external meshing pair at different speeds.

Figure 4 Vibration response of the internal meshing pair at different speeds.

The time history diagram and phase plane diagram of the external meshing pair at n=2000, n=4000 and n=8000 rpm (revolutions per minute) is shown in Fig. 3. The vibration response of the system at each speed shows the period responses of 3 times, 2 times and 1 time, respectively. Similarly, at three input speeds, the time history diagram and phase plane diagram of the internal meshing pair in Fig. 4 show similar vibration response laws, and the vibration response tends, from multiple periods, to form a single period vibration.

Figure 5 shows the vibration displacement frequency response curve of the external and internal meshing pairs with the change in speed. The vertical axis, X_max and X_min, are the maximum and minimum vibration displacement, respectively, and the horizontal axis, n, is the input speed.

Figure 5 Vibration response curve of the planetary gear system.

Figure 5a shows that under the action of various kinds of excitation, including clearance, the vibration displacement frequency response of the external meshing pair jumps at n=9100 rpm, the system transits from a non-impact state to unilateral impact state, and the system is in unilateral impact state in the speed range of n = 9100–9700 rpm. Figure 5b shows that the vibration displacement frequency response of the internal meshing pair jumps at n=7100 rpm, and the range of unilateral impact speed is n = 7100–7700 rpm. At the same time, in addition to the main resonance peak, the frequency response curves of the internal and external meshing pairs also form the resonance sub-peak, which leads to the super harmonic vibration response of the system.

4.2 Influence of excitation parameters on frequency response characteristics

4.2.1 Meshing stiffness influence

Based on the basic and dynamic parameters in Tables 1 and 2, taking ${\mathit{\xi }}_{\mathrm{sp}\mathrm{1}}={\mathit{\xi }}_{\mathrm{rp}\mathrm{1}}=\mathrm{0.035}$, the other excitation parameters are L=3, N=3, F₁=1, f₁=0.05, f₂=0.02, f₃=0.01, e₁=0.05, e₂=0.02, and e₃=0.01. The vibration displacement response curves under the three kinds of time-varying meshing stiffness are obtained and shown in Fig. 6. The harmonic components of the time-varying meshing stiffness are k₁=0, k₂=0, k₃=0, k₁=0.05, k₂=0.02, k₃=0.01 and k₁=0.1, k₂=0.04, and k₃=0.05.

Figure 6 Effect of stiffness on the vibration response curve.

Figure 6 shows that the harmonic component of the variable meshing stiffness decreases, the jumping phenomenon of the vibration response curves of the external and internal meshing pairs decreases, the speed range of the system with unilateral impact decreases, and the nonlinear vibration of the system is significantly suppressed.

4.2.2 Meshing damping influence

Given that the meshing damping ratio ξ_sp1 is equal to 0.025, 0.04 and 0.1, respectively, we define ξ_rp1=ξ_sp1. In addition, the harmonic component of external meshing is k₁=0.1, k₂=0.04, and k₃=0.02. The internal harmonic component is k₁=0.05, k₂=0.02, and k₃=0.01. The other excitation parameters are L=3, N=3, F₁=1, f₁=0.05, f₂=0.02, f₃=0.01, e₁=0.05, e₂=0.02, and e₃=0.01. Based on the above parameters, the solution result of the vibration displacement response curve of the system is shown in Fig. 7.

Figure 7 Effect of damping on the vibration response curve.

Figure 7 shows that the nonlinear jump phenomenon in the vibration displacement response curves of the external and internal meshing pairs decreases with the increase in meshing damping and, finally, disappears. There is no unilateral impact, and the system vibration presents linear characteristics.

4.2.3 External exciting force influence

When the external exciting force F₁ is equal to 1, 1.25 and 1.5, respectively, the vibration displacement response curve of the external and internal meshing pairs is as shown in Fig. 8. The other parameters are L=3, N=3, f₁=0.05, f₂=0.02, f₃=0.01, ξ_sp1=0.025, ξ_rp1=0.025, e₁=0.05, e₂=0.02, e₃=0.01, k₁=0.1, k₂=0.04, and k₃=0.02.

Figure 8 Effect of the external exciting force on the vibration response curve.

Figure 8 shows that the amplitude of the vibration displacement response curve of the external and internal meshing pairs increases obviously with the increase in the external exciting force, while the nonlinear jump phenomenon and unilateral impact area remain basically unchanged, and the nonlinear characteristics of the system are less affected by an external exciting force.

4.2.4 Backlash influence

Figure 9 shows the speed vibration displacement response curves of the internal and external meshing pairs under different backlash and taking $b_{\mathrm{sp}\mathrm{1}}=b_{\mathrm{rp}\mathrm{1}}=\mathrm{10}$ µm, $b_{\mathrm{sp}\mathrm{1}}=b_{\mathrm{rp}\mathrm{1}}=\mathrm{5}$ µm, and $b_{\mathrm{sp}\mathrm{1}}=b_{\mathrm{rp}\mathrm{1}}=\mathrm{0}$ µm. The other excitation parameters can be expressed as L=3, N=3, F₁=1, f₁=0.05, f₂=0.02, f₃=0.01, ξ_sp1=0.025, ξ_rp1=0.025, e₁=0.05, e₂=0.02, e₃=0.01, k₁=0.1, k₂=0.04, and k₃=0.02.

Figure 9 shows that when the backlash is zero, the vibration response curve of the system is linear. When the backlash is not zero, the vibration displacement response curves of the external and internal meshing pairs appear to have a jumping phenomenon and a unilateral impact area. At the same time, with the increase in the backlash, the nonlinear characteristics of the system are enhanced, and compared with zero backlash, the main resonance peak of the vibration response curve of the system is suppressed.

Figure 9 Effect of backlash on the vibration response curve.

5 Frequency response vibration test

The vibration test bench is shown in Fig. 10. The main components include the electric control cabinet, data acquisition instrument, drive motor, planetary gearbox, acceleration sensor, torque sensor and magnetic particle brake. The basic parameters of the planetary gearbox are shown in Table 1.

Figure 10 Vibration test bench.

The characteristic frequency is calculated with Eqs. (9), (10) and (11) and is shown in Table 3.

Figures 11, 12 and 13 are the time domain and frequency domain diagrams of the vibration acceleration signal of the planetary gearbox when the motor selects the uniform speed mode of 1500, 2000 and 2500 rpm, respectively.

Figure 11 Vibration acceleration signal of the gearbox at n=1500 rpm.

Figure 12 Vibration acceleration signal of the gearbox at n=2000 rpm.

Table 3 Characteristic frequency of the planetary gearbox.

Figure 13 Vibration acceleration signal of the gearbox at n=2500 rpm.

Figures 11b, 12b and 13b show that, under three input speeds, the frequency domain signal of vibration acceleration of planetary gearbox has obvious vibration peaks at three characteristic frequencies of the gearbox, which indicates that the vibration acceleration signal measured in the experiment has a good level of reliability.

The time domain signals of the vibration displacement of the planetary gearbox can be obtained by processing the time domain signals of the three groups of the gearbox vibration acceleration signals measured by the analysis software, as shown in Fig. 14. At the same time, based on the system parameters of the planetary gearbox, the theoretical spectrum diagram of the system under three input speeds is obtained through theoretical analysis.

Figure 14 Time domain vibration displacement signals of the gearbox at different speeds.

Figure 15 shows the spectrum of the theoretical and experimental displacement of the system when the input speed n=1500 rpm. The calculated meshing frequency of the system is 360 Hz at this speed. Figure 15a is the result of theoretical analysis; the black solid line is the spectrum curve of the external meshing pair under the input speed in the theoretical analysis, and the red solid line is the spectrum curve of the internal meshing pair. The analysis of the test results is shown in Fig. 15b.

Compared with Fig. 15a and b, the results show that, in the theoretical analysis, both the internal and external meshing pairs of the system produce four resonance peaks, and the peak frequency points are the same, which are, respectively, close to the system one-quarter frequency multiplication, one-half frequency multiplication, three-quarters frequency multiplication and the meshing frequency points. In addition to the above four points, the spectrum curve of the system is smooth without side frequency resonance. In the experimental analysis, the planetary gearbox also resonates near the frequency point of the period doubling, and there is also side frequency resonance near the frequency point. In addition to the resonance caused by the motor input speed, there are other forms of vibration signals in the system.

Figure 15 Frequency domain diagram of the vibration displacement at n=1500 rpm.

Figure 16 Frequency domain diagram of the vibration displacement at n=2000 rpm.

Figure 17 Frequency domain diagram of the vibration displacement at n=2500 rpm.

When n=1500 rpm, the corresponding frequency values of frequency multiplication response of the planetary gearbox in the theoretical analysis and experimental analysis are shown in Table 4.

Table 4 Frequency multiplication response of the gearbox at n=1500 rpm.

Table 5 Frequency multiplication response of the gearbox at n=2000 rpm.

Table 6 Frequency multiplication response of the gearbox at n=2500 rpm.

Figure 16 shows the spectrum of the theoretical and experimental displacement of the system at the input speed of n=2000 rpm. The calculated meshing frequency of the system at this speed is 480 Hz. Compared with Fig. 16a and b, the results show that the theoretical analysis spectrum produces a resonance response at four frequency points during the frequency multiplication period, while the experimental analysis spectrum only forms resonance peak at the one-quarter frequency multiplication point and one-half frequency multiplication point.

When n=2000 rpm, the corresponding frequency values of the frequency multiplication response of the planetary gearbox in theoretical and experimental analysis are shown in Table 5.

Figure 17 shows the spectrum of the theoretical and experimental displacement of the system at an input speed of n=2500 rpm. The calculated meshing frequency of the system is 600.12 Hz at this speed. Compared with Fig. 17a and b, the results show that the theoretical spectrum produces a resonance response at four frequency points during the frequency multiplication period, while the experimental spectrum only forms a resonance peak at the one-quarter frequency multiplication point, one-half frequency multiplication point and the meshing frequency.

When n=2500 rpm, the corresponding frequency values of the frequency multiplication response of the planetary gearbox in the theoretical analysis and experimental analysis are shown in Table 6.

By analyzing the spectrum diagrams under three input speeds, it can be determined that the amplitude of the test spectrum diagram is larger than that of the theoretical spectrum diagram, and the vibration signal of the test gearbox itself is more complex. Therefore, the difference in amplitude can be regarded as reasonable. At the same time, the response points of each frequency multiple of the theoretical and experimental spectrum diagram are similar, and the characteristics of vibration spectrum are the same and have similar vibration response laws. Therefore, the correctness of the analytical model established by the lumped parameter method and the effectiveness of the theoretical analysis results can be confirmed.

6 Conclusion

In order to solve the problem of the nonlinear frequency response of planetary gear transmission system, based on the lumped parameter theory, a nonlinear frequency response model of the bending torsion coupling of a planetary gear system is established by comprehensively considering the backlash, support clearance, time-varying meshing stiffness, transmission error, external periodic excitation and other factors. The Runge–Kutta numerical method is used to solve the model, and the nonlinear frequency response characteristics of the system are analyzed. Finally, the vibration test of the planetary gearbox is carried out. The conclusions are as follows:

1. With the increase in input speed, the vibration response of the planetary gear system transits from a multiple-period response to a single-period response. Under the joint action of backlash, time-varying meshing stiffness and meshing damping, the speed vibration displacement response curve of the internal and external meshing pair appears to have a nonlinear jump phenomenon and unilateral impact area, and the system presents strong nonlinear characteristics.

2. The nonlinear vibration can be effectively suppressed by decreasing the time-varying meshing stiffness or increasing the meshing damping. However, by increasing the external excitation force, the vibration response displacement of the system increases, and the nonlinear characteristics of the system remain basically unchanged. The existence of backlash is an important inducement to the nonlinear vibration of the system, but it can inhibit the vibration of the system in a certain range.

3. In the vibration test, the vibration spectrum characteristics of the vibration displacement frequency domain signal of the planetary gearbox at different speeds are similar to the simulation results, which proves the validity of the simulation analysis model and simulation results.