Semi-numerical analysis of a two-stage series composite planetary transmission considering IHB and MsP methods

Wang, Yongmei; Wang, Xigui; Fu, Baixue; Ruan, Jiafu; An, Siyuan

doi:10.21595/jve.2021.21832

Journal of Vibroengineering

Browse Journal

Submit article

Published: 06 May 2021

Check for updates

Semi-numerical analysis of a two-stage series composite planetary transmission considering IHB and MsP methods

Yongmei Wang¹

Xigui Wang²

Baixue Fu³

Jiafu Ruan⁴

Siyuan An⁵

^{1, 3}School of Motorcar Engineering, Heilongjiang Institute of Technology, No. 999, Hongqidajie Road, Daowai District, Harbin, 150036, P. R. China

^{2, 4, 5}School of Engineering Technology, Northeast Forestry University, No. 26, Hexing Road, Xiangfang District, Harbin, 150040, P. R. China

Corresponding Author:

Xigui Wang

Cite the article Download PDF

Downloads 829

Abstract

This study conducts an analytical investigation of the dynamic response characteristics of a two-stage series composite system (TsSCS) with a planetary transmission consisting of dual-power branches. An improved incremental harmonic balance (IHB) method, based on the arc length extension technique, is proposed. The results are compared with those of the numerical integration method to verify the feasibility and effectiveness of the improved method. Following that, the multi-scale perturbation (MsP) method is applied to the TsSCS proposed in this paper. The frequency response equations of the primary resonance, subharmonic resonance, and superharmonic resonance are solved to generate the frequency response characteristic curves of the planetary transmission system in this method. A comparison between the results obtained by the MsP method and the numerical integration method proves that the former is ideal and credible in most regions. Considering the parameters of TsSCS excitation frequency and damping, the nonlinear response characteristics of steady-state motion are mutually converted. The effects of the time-varying parameters and the nonlinear deenthing caused by the gear teeth clearance on the amplitude-frequency characteristics of TsSCS components are studied in this special topic.

Semi-numerical analysis of a two-stage series composite planetary transmission considering IHB and MsP methods

Highlights

The dynamic response characteristics of a two-stage series composite system with a planetary transmission system that consists of dual-power branches are analyzed.
An alternative to the computationally heavy and financially expensive numerical methods to analyze planetary transmission systems is proposed.
A comparison between the results by the MsP method and the numerical integration method proves that the former is ideal and credible in most regions are obtained.
The parameters of TsSCS excitation frequency and damping, the nonlinear response characteristics of steady-state motion are mutually converted.
The effects of the time-varying parameters and the nonlinear deenthing caused by the gear teeth clearance on the amplitude-frequency characteristics of TsSCS components are studied.

1. Introduction

Planetary gears are the most critical underwater components of ships that transfer real-time power and time-varying motion. Owing to its compact transmission structure, strong anti-scuffing bearing performance, and high transmission precision, the planetary gear transmission is widely used in various mechanical systems. The application of the planetary transmission system in the underwater military industry has created an urgent demand for lightweight and reliable structures [1]. Planetary gears have various dynamic response properties due to their complex inherent characteristics [2]. The dynamic behavior analysis of the alternating meshing process is a popular area of study in the design and applications of various planetary gear mechanisms and their transmission systems [3-5]. Many researchers have been working on new simulation analysis methods to gain an in-depth understanding of the dynamic meshing transients of planetary gears [6-8]. Numerous studies have been performed on the dynamic behavior analysis method, which plays an important role in determining the performance of planetary gears. In addition, the machinery industry is paying more attention to the dynamic characteristics of planetary gear drives and the resulting vibrations and noise [9-11].

A two-stage series composite system (TsSCS), consisting of a planetary transmission system with dual-power branches, is a complex and flexible mechanical system that comprises several parts [12]. A TsSCS is generally divided into two parts: the transmission system (gear train, transmission shaft, and bearing) and the structural system (gearbox, dual-branch composite planetary gearbox, dynamometer, and bracket) [13]. The increasing number of studies on the vibrations and noise produced by underwater devices, particularly on controlling the noise produced by the power rear transmission system of underwater devices, must account for the TsSCS and its meshing process, in addition to the real-time dynamic meshing force acting on the system. Fig. 1 reveals outline drawing and cutaway view of a certain type of warship TsSCS.

Fig. 1Outline drawing and cutaway view of a certain type of TsSCS

a) Outline drawing of a certain type of TsSCS

b) Cutaway view of a certain type of TsSCS

Few theoretical studies consider the TsSCS with a planetary transmission consisting of dual-power branches as a subject [14]. It is important to possess some professional design knowledge while studying the unique kinematics and geometric characteristics of a TsSCS consisting of a planetary transmission with dual-power branches [15]. The TsSCS, with a dual-power branch planetary transmission, is preferred to the horizontal shaft gear deceleration system, especially in applications requiring high linear speed power density designs and kinematic flexibility to optimize different speed ratios [16-18]. It has been demonstrated that reducing the spoke thickness to increase gear flexibility also resolves several internal gear and planetary frame errors and operational errors, in addition to making the system lighter [19]. Moreover, a flexible internal gear improves the load sharing between planets, which is an important feature if manufacturing and assembly related gear and carrier errors are inevitable. Thus, it is difficult to quantify the factors that influence the TsSCS with a dual-power branch planetary transmission under quasi-static conditions. The research in this subject is TsSCS with a planetary transmission consisting of dual-power branches, which uses a translation-torsion coupling model. It is difficult to solve the problem of dynamic differential equations by applying existing analytical methods. To deeply study the dynamic characteristics of TsSCS, it is necessary to find a suitable new method to solve this kind of complex model. This article mainly studies the analytical solution method of TsSCS kinetic equation. The translation-torsion coupling dynamic equation of TsSCS with a planetary transmission consisting of dual-power branches is non-dimensionalized.

Modals are the inherent characteristics of gear transmission systems [20-22]. A modal analysis is used to determine the vibration characteristics of the designed structure or its transmission components. The modal analysis is a part of the structural dynamics analysis and is also the starting point for the subsequent transient dynamics, harmonic response, and spectrum analyses [23]. Each mode has a corresponding natural frequency, damping ratio, and mode shape [24]. A modal analysis is a modern technique that is used to study the dynamic characteristics of a transmission system structure, including modal analysis of linear vibration theory and experimental modal analysis [25]. An experimental modal analysis can only be carried out after the components of the structure have been processed and assembled. However, it has a longer test cycle and is more expensive. In addition, it is easily affected by the quality of the processing and assembly steps. Thus, it is difficult to use these results in the design analysis stage [26]. As a result, the experimental analysis is used to verify the results of the analysis of the theoretical model and modify it accordingly. Various modal parameters, such as the modal frequency, shape, quality, stiffness, damping, etc., affect the dynamic load design.

Since most of the internal, planetary, and sun gears are excluded from these models, it is not possible to study the effect of the internal gear thickness on the performance of the TsSCS and its effect on the stress acting on the planetary and sun gears and the load sharing between the planets. Similarly, it is not possible to accurately predict the shape and deflection of the gears. Although the above mentioned studies indicate that the adverse effects of the gear and planet carrier manufacturing errors can be minimized by improving the planetary load-sharing characteristics under quasi-static conditions, these modifications lead to increased gear contact stress. These static analyses alone cannot predict the actual design of the system because the increasing flexibility of the TsSCS causes its performance to change only under dynamic conditions. This might also lead to an rise in the stress acting on the gear. The current study, thus, studies the inherent characteristics of the TsSCS, constructs its dynamic model, and synthetically analyzes its inherent and dynamic response characteristics.

2. Mathematical model of the modal analysis

The first step of the dynamic calculation is to determine the natural frequency and mode shape of the structure while ignoring its damping. These results reflect the basic dynamic characteristics of the structure and its response trend under dynamic loads. The double-branch compound transmission system test bench is considered to be a single elastic body. The differential equation of motion describing the overall vibration of the test bench is given as:

1

M \ddot{u} (t) + C \dot{u} (t) + K u (t) = Q (t),

where $\ddot{u} (t)$ , $\dot{u} (t)$ , and $u (t)$ refer to the acceleration, velocity, and displacement vectors of the nodes in the vibration system, respectively. $M$ , $C$ , and $K$ represent the mass, damping, and stiffness matrices of the vibration system, respectively. $Q (t)$ represents the external force vector received by node. Only the inherent characteristics of the vibration system are modeled and solved in the modal analysis. The model does not contain external force terms and neglects the damping terms that are assumed to have an insignificant impact on the system. The differential equation of motion for an undamped free vibration system is given as:

2

M \ddot{u} (t) + K u (t) = 0 .

The resonance solution form of the equation is given as:

3

u (t) = ϕ s i n ω_{n} t,

where $u$ is the displacement vector and $ϕ$ is the characteristic vector that amplitude of the displacement vector $u$ . $ω_{n}$ is the natural angular frequency.

The resonance form of the equation is a key to the numerical solution. It is derived under the assumption that all degrees of freedom of the vibrating structure move in a synchronous manner. During this process, the basic shape of the structure does not change, only the amplitude varies.

2.1. Dimensionless differential equation of the second-stage star gear train

(1) Dimensionless differential equation of the ring gear of the star gear train:

4

\{\begin{array}{l} m_{r}^{Ι} {\hat{ω}}^{2} {\ddot{y}}_{1} - \sum_{m = 1}^{4} (\hat{ω} C_{r m}^{Ι} {\dot{\hat{δ}}}_{r m}^{Ι} + k_{r m}^{Ι} {\hat{δ}}_{r m}^{Ι}) s i n ψ_{r m} c o s β_{b m} + (\hat{ω} C_{r}^{Ι} {\dot{y}}_{1} + k_{r}^{Ι} y_{1}) = 0, \\ m_{r}^{Ι} {\hat{ω}}^{2} {\ddot{y}}_{2} + \sum_{m = 1}^{4} (\hat{ω} C_{r m}^{Ι} {\dot{\hat{δ}}}_{r m}^{Ι} + k_{r m}^{Ι} {\hat{δ}}_{r m}^{Ι}) c o s ψ_{r m} c o s β_{b m} + (\hat{ω} C_{r}^{Ι} {\dot{y}}_{2} + k_{r}^{Ι} y_{2}) = 0, \\ (I_{p}^{Ι Ι} / (r_{p}^{Ι Ι})^{2}) {\hat{ω}}^{2} {\ddot{y}}_{3} + \sum_{m = 1}^{4} (\hat{ω} C_{r m}^{Ι} {\dot{\hat{δ}}}_{r m}^{Ι} + k_{r m}^{Ι} {\hat{δ}}_{r m}^{Ι}) c o s β_{b m} + (\hat{ω} C_{r u}^{Ι} {\dot{y}}_{3} + k_{r u}^{Ι} y_{3}) \\ + \frac{k_{r s u}^{12}}{r_{r}^{Ι}} (\frac{y_{3}}{r_{r}^{Ι}} - \frac{y_{24}}{r_{s}^{Ι Ι}}) + \frac{\hat{ω} C_{r s u}^{12}}{r_{r}^{Ι}} (\frac{{\dot{y}}_{3}}{r_{r}^{Ι}} - \frac{{\dot{y}}_{24}}{r_{s}^{Ι Ι}}) = 0 . \end{array}

(2) Dimensionless differential equation of the sun gear of the star gear train:

5

\{\begin{array}{l} m_{s}^{Ι} {\hat{ω}}^{2} {\ddot{y}}_{4} - \sum_{m = 1}^{4} (\hat{ω} C_{s m}^{Ι} {\dot{\hat{δ}}}_{s m}^{Ι} + k_{s m}^{Ι} {\hat{δ}}_{s m}^{Ι}) s i n ψ_{s m} c o s β_{b m} + (\hat{ω} C_{s}^{Ι} {\dot{y}}_{4} + k_{s}^{Ι} y_{4}) = 0, \\ m_{s}^{Ι} {\hat{ω}}^{2} {\ddot{y}}_{5} + \sum_{m = 1}^{4} (\hat{ω} C_{s m}^{Ι} {\dot{\hat{δ}}}_{s m}^{Ι} + k_{s m}^{Ι} {\hat{δ}}_{s m}^{Ι}) c o s ψ_{s m} c o s β_{b m} + (\hat{ω} C_{s}^{Ι} {\dot{y}}_{5} + k_{s}^{Ι} y_{5}) = 0, \\ (I_{p}^{Ι Ι} / (r_{p}^{Ι Ι})^{2}) {\hat{ω}}^{2} {\ddot{y}}_{6} + \sum_{m = 1}^{4} (\hat{ω} C_{s m}^{Ι} {\dot{\hat{δ}}}_{s m}^{Ι} + k_{s m}^{Ι} {\hat{δ}}_{s m}^{Ι}) c o s β_{b m} + (\hat{ω} C_{s u}^{Ι} {\dot{y}}_{6} + k_{s u}^{Ι} y_{6}) = \frac{T_{s}^{Ι}}{r_{s}^{Ι} {\bar{e}}_{s 1}^{Ι}} . \end{array}

(3) Dimensionless differential equation of the star gear of the star gear train:

6

\{\begin{array}{l} m_{p}^{Ι} {\hat{ω}}^{2} {\ddot{y}}_{7 + 3 (m - 1)} - (\hat{ω} C_{s m}^{Ι} {\dot{\hat{δ}}}_{s m}^{Ι} + k_{s m}^{Ι} {\hat{δ}}_{s m}^{Ι}) s i n α_{s m} c o s β_{b m} \\ + (\hat{ω} C_{r m}^{Ι} {\dot{\hat{δ}}}_{r m}^{Ι} + k_{r m}^{Ι} {\hat{δ}}_{r m}^{Ι}) s i n α_{r m} c o s β_{b m} + (\hat{ω} C_{p m}^{Ι} {\dot{y}}_{7 + 3 (m - 1)} + k_{p m}^{Ι} y_{7 + 3 (m - 1)}) = 0, \\ m_{p}^{Ι} {\hat{ω}}^{2} {\ddot{y}}_{8 + 3 (m - 1)} - (\hat{ω} C_{s m}^{Ι} {\dot{\hat{δ}}}_{s m}^{Ι} + k_{s m}^{Ι} {\hat{δ}}_{s m}^{Ι}) c o s α_{s m} c o s β_{b m} \\ - (\hat{ω} C_{r m}^{Ι} {\dot{\hat{δ}}}_{r m}^{Ι} + k_{r m}^{Ι} {\hat{δ}}_{r m}^{Ι}) c o s α_{r m} c o s β_{b m} + (\hat{ω} C_{p m}^{Ι} {\dot{y}}_{8 + 3 (m - 1)} + k_{p m}^{Ι} y_{8 + 3 (m - 1)}) = 0, \\ (I_{p}^{Ι} / (r_{p}^{Ι})^{2}) {\hat{ω}}^{2} {\ddot{y}}_{9 + 3 (m - 1)} + (\hat{ω} C_{s m}^{Ι} {\dot{\hat{δ}}}_{s m}^{Ι} + k_{s m}^{Ι} {\hat{δ}}_{s m}^{Ι}) c o s β_{b m} \\ - (\hat{ω} C_{r m}^{Ι} {\dot{\hat{δ}}}_{r m}^{Ι} + k_{r m}^{Ι} {\hat{δ}}_{r m}^{Ι}) c o s β_{b m} + (\hat{ω} C_{p m u}^{Ι} {\dot{y}}_{9 + 3 (m - 1)} + k_{p m u}^{Ι} y_{9 + 3 (m - 1)}) = 0 . \end{array}

On the basis of the non-dimensionalized equation, the incremental harmonic balance method and the multi-scale perturbation analysis method based on the arc-length continuation technology are respectively used to obtain It is suitable for the solution of TsSCS kinetic equation.

2.2. Dimensionless differential equation of the second-stage planetary gear train

(1) Dimensionless differential equation of the planet carrier of the planetary gear train:

7

\{\begin{array}{l} m_{c}^{Ι Ι} {\hat{ω}}^{2} {\ddot{y}}_{19} + \sum_{n = 1}^{4} [k_{p n}^{Ι Ι} ({\hat{δ}}_{n r}^{Ι Ι} c o s ψ_{n} - {\hat{δ}}_{n t}^{Ι Ι} s i n ψ_{n}) + \hat{ω} C_{p n}^{Ι Ι} ({\dot{\hat{δ}}}_{n r}^{Ι Ι} c o s ψ_{n} - {\dot{\hat{δ}}}_{n t}^{Ι Ι} s i n ψ_{n})] \\ + (\hat{ω} C_{c}^{Ι Ι} {\dot{y}}_{19} + k_{c}^{Ι Ι} y_{19}) = 0, \\ m_{c}^{Ι Ι} {\hat{ω}}^{2} {\ddot{y}}_{20} + \sum_{n = 1}^{4} [k_{p n}^{Ι Ι} ({\hat{δ}}_{n r}^{Ι Ι} s i n ψ_{n} + {\hat{δ}}_{n t}^{Ι Ι} c o s ψ_{n}) + \hat{ω} C_{p n}^{Ι Ι} ({\dot{\hat{δ}}}_{n r}^{Ι Ι} s i n ψ_{n} + {\dot{\hat{δ}}}_{n t}^{Ι Ι} c o s ψ_{n})] \\ + (\hat{ω} C_{c}^{Ι Ι} {\dot{y}}_{20} + k_{c}^{Ι Ι} y_{20}) = 0, \\ (I_{c}^{Ι Ι} / r_{c}^{Ι Ι})^{2} {\hat{ω}}^{2} {\ddot{y}}_{21} + \sum (k_{p n}^{Ι Ι} {\hat{δ}}_{n t}^{Ι Ι} + \hat{ω} C_{p n}^{Ι Ι} {\dot{\hat{δ}}}_{n t}^{Ι Ι}) + (\hat{ω} C_{c u}^{Ι Ι} {\dot{y}}_{21} + k_{c u}^{Ι Ι} y_{21}) = T_{c}^{Ι Ι} / r_{c}^{Ι Ι} {\bar{e}}_{s 1}^{Ι} . \end{array}

(2) Dimensionless differential equation of the sun gear of the planetary gear train:

8

\{\begin{array}{l} m_{s}^{Ι Ι} {\hat{ω}}^{2} {\ddot{y}}_{22} - \sum_{n = 1}^{4} (\hat{ω} C_{s n}^{Ι Ι} {\dot{\hat{δ}}}_{s n}^{Ι Ι} + k_{s n}^{Ι Ι} {\hat{δ}}_{s n}^{Ι Ι}) s i n ψ_{s n} c o s β_{b n} + (\hat{ω} C_{s}^{Ι Ι} {\dot{y}}_{22} + k_{s}^{Ι Ι} y_{22}) = 0, \\ m_{s}^{Ι Ι} {\hat{ω}}^{2} {\ddot{y}}_{23} + \sum_{n = 1}^{4} (\hat{ω} C_{s n}^{Ι Ι} {\dot{\hat{δ}}}_{s n}^{Ι Ι} + k_{s n}^{Ι Ι} {\hat{δ}}_{s n}^{Ι Ι}) c o s ψ_{s n} c o s β_{b n} + (\hat{ω} C_{s}^{Ι Ι} {\dot{y}}_{23} + k_{s}^{Ι Ι} y_{23}) = 0, \\ (I_{s}^{Ι Ι} / (r_{s}^{Ι Ι})^{2}) {\hat{ω}}^{2} {\ddot{y}}_{24} + \sum_{n = 1}^{4} (\hat{ω} C_{s n}^{Ι Ι} {\dot{\hat{δ}}}_{s n}^{Ι Ι} + k_{s n}^{Ι Ι} {\hat{δ}}_{s n}^{Ι Ι}) c o s β_{b n} + (\hat{ω} C_{s u}^{Ι Ι} {\dot{y}}_{24} + k_{s u}^{Ι Ι} y_{24}) \\ + [\hat{ω} C_{r s u}^{12} ({\dot{y}}_{24} / r_{s}^{Ι Ι} - {\dot{y}}_{3} / r_{r}^{Ι}) + k_{r s u}^{12} (y_{24} / r_{s}^{Ι Ι} - y_{3} / r_{r}^{Ι})] / r_{s}^{Ι Ι} = T_{s}^{Ι Ι} / r_{s}^{Ι Ι} {\bar{e}}_{s 1}^{Ι} \end{array}

(3) Dimensionless differential equation of the planetary gear of the planetary gear train:

9

\{\begin{array}{l} m_{p}^{Ι Ι} {\hat{ω}}^{2} {\ddot{y}}_{25 + 3 (n - 1)} - (\hat{ω} C_{s n}^{Ι Ι} {\dot{\hat{δ}}}_{s n}^{Ι Ι} + k_{s n}^{Ι Ι} {\hat{δ}}_{s n}^{Ι Ι}) s i n α_{s n} c o s β_{b n} \\ + (\hat{ω} C_{r n}^{Ι Ι} {\dot{\hat{δ}}}_{r n}^{Ι Ι} + k_{r n}^{Ι Ι} {\hat{δ}}_{r n}^{Ι Ι}) s i n α_{r n} c o s β_{b n} - (\hat{ω} C_{p n}^{Ι Ι} {\dot{\hat{δ}}}_{n r}^{Ι Ι} + k_{p n}^{Ι Ι} {\hat{δ}}_{n r}^{Ι Ι}) = 0, \\ m_{p}^{Ι Ι} {\hat{ω}}^{2} {\ddot{y}}_{26 + 3 (n - 1)} - (\hat{ω} C_{s n}^{Ι Ι} {\dot{\hat{δ}}}_{s n}^{Ι Ι} + k_{s n}^{Ι Ι} {\hat{δ}}_{s n}^{Ι Ι}) c o s α_{s n} c o s β_{b n} \\ - (\hat{ω} C_{r n}^{Ι Ι} {\dot{\hat{δ}}}_{r n}^{Ι Ι} + k_{r n}^{Ι Ι} {\hat{δ}}_{r n}^{Ι Ι}) c o s α_{r n} c o s β_{b n} - (\hat{ω} C_{p n}^{Ι Ι} {\dot{\hat{δ}}}_{n r}^{Ι Ι} + k_{p n}^{Ι Ι} {\hat{δ}}_{n t}^{Ι Ι}) = 0, \\ (I_{p}^{Ι Ι} / (r_{p}^{Ι Ι})^{2}) {\hat{ω}}^{2} {\ddot{y}}_{27 + 3 (n - 1)} + (\hat{ω} C_{s n}^{Ι Ι} {\dot{\hat{δ}}}_{s n}^{Ι Ι} + k_{s n}^{Ι Ι} {\hat{δ}}_{s n}^{Ι Ι}) c o s β_{b n} - (\hat{ω} C_{r n}^{Ι Ι} {\dot{\hat{δ}}}_{r n}^{Ι Ι} + k_{r n}^{Ι Ι} {\hat{δ}}_{r n}^{Ι Ι}) c o s β_{b n} \\ - (\hat{ω} C_{p u n}^{Ι Ι} {\dot{y}}_{27 + 3 (n - 1)} + k_{p u n}^{Ι Ι} y_{27 + 3 (n - 1)}) = 0 . \end{array}

Rewriting Eq. (9) in its matrix form:

10

\hat{M} \ddot{y} + \hat{c} (\dot{y}, τ) + \hat{h} (y, τ) = \hat{M} \ddot{y} + \hat{C} (\dot{y}, τ) \dot{y} + \hat{K} (y, ω_{m i} τ) y = \hat{F} (τ) .

It is noteworthy that the stiffness matrix $\hat{K}$ is no longer symmetric due to the transient nature of the phase angle of the planetary gears.

3. Incremental harmonic balance and arc length extension method

The harmonic balance method uses a description function to approximate the nonlinearity caused by the gap and is widely used in planetary transmission systems. The excitation and response parameters are assumed to be harmonic functions and are substituted in the nonlinear equation. The approximate expressions of the response and phase parameters can thus be obtained by using the condition of equal power coefficients. Since this method is not limited by the degree of nonlinearity, all the frequency response values can be obtained. However, owing to the limitations of the assumed excitation and response form, the accuracy of this method is not satisfactory, particularly if the first harmonic is considered. It results in the artificial loss of the superharmonic, subharmonic, or chaotic responses. Lau and Cheung proposed the incremental harmonic balance (IHB) method in 1981 to improve the accuracy of the existing harmonic balance method. A Taylor series expansion was performed on the nonlinear differential equations while ignoring the higher-order derivatives to obtain the differential equations in an incremental form. The Fourier series and the Galerkin method were then used to obtain nonlinear algebraic equations. The entire process is divided into an incremental component (Newton-Raphson method) and a harmonic balance component (Galerkin method). This method has the advantage of free control algorithm convergence accuracy, among many others, and is thus an effective method for solving complex nonlinear problems [27]. Rohan and Lukeš used the incremental harmonic balance method for a two-stage star gear train with multiple degrees of freedom; however, it is only used the response as an incremental parameter [28]. If a singular point is encountered along the path of the solution branch, the quasi-arc length parameter is introduced. The original variables and parameters are assumed to be functions of the arc length. This condition is added to the original equation to smoothly track the path through the singular point [29]. The harmonic balance method, based on the continuation of the arc-length, has been applied to the pure torsion model of a fixed shaft and a single-stage planetary transmission [30]. The incremental harmonic balance method, based on the continuation of the arc length, is used to calculate the dynamic characteristic equation of the system used in this study. The two incremental parameters (response and fundamental frequency) are expressed in terms of the arc length to smoothly overcome the singular points along the path. The formula to calculate the steady-state response of a two-stage herringbone planetary transmission system is presented in this study. This method has not yet been applied to planetary gear transmission systems.

3.1. Incremental harmonic balance method

A new time variable $τ_{h} = Ω τ$ is introduced in this method. The expression of $y$ , initially written in terms of $τ$ , is rewritten in terms of $τ_{h}$ . As a result, Eq. (10) is written as:

11

Ω^{2} \hat{M} y^{″} + Ω \hat{C} (\dot{y}, τ_{h}) y^{'} + \hat{K} (y, \frac{ω_{m i}}{Ω} τ_{h}) y = \hat{F} (τ_{h}) .

The incremental process is the first component of the IHB method. If $y_{j 0}$ and $Ω_{0}$ are the solutions of Eq. (11), their neighboring points can be expressed as:

12

y_{j} = y_{j 0} + Δ y_{j}, Ω = Ω_{0} + Δ Ω, ​ ​ ​ ​ ​ ​ ​ ​ ​ ​ ​ ​ ​ ​ ​ ​ ​ ​ ​ j = 1,2, \dots, M + N,

where $Δ y_{j}$ and $Δ Ω$ are the incremental parameters.

By substituting Eq. (12) into Eq. (11) and omitting high-order small quantities, the incremental equation matrix can be obtained with $Δ y_{j}$ and $Δ Ω$ as the unknown quantities:

13

Ω_{0}^{2} \hat{M} Δ y^{″} + Ω_{0} \hat{C} Δ y^{'} + \hat{K} Δ y = R - (2 Ω_{0} \hat{M} y_{0}^{″} + \hat{C} y_{0}^{'}) Δ Ω,

14

R = \hat{F} (τ_{h}) - [Ω_{0}^{2} \hat{M} y_{0}^{″} + Ω_{0} \hat{C} y_{0}^{'} + \hat{K} y_{0}],

where $R$ is the unbalanced force vector (also referred to as the residual correction term in some studies). If $y_{j 0}$ and $Ω_{0}$ are exact solutions, then $R = 0$ .

The second component of the IHB method involves the harmonic balance process. The steady-state response of the system is described by a Fourier series. The response contains only odd harmonics and is given as follows:

15

y_{j 0} = a_{j 0} + \sum_{k = 1}^{N_{c}} [a_{j k} \cos k τ_{h} + b_{j k} \sin k τ_{h}] = C_{s} A_{j},

16

Δ y_{j 0} = Δ a_{j 0} + \sum_{k = 1}^{N_{c}} [Δ a_{j k} \cos k τ_{h} + Δ b_{j k} \sin k τ_{h}] = C_{s} Δ A_{j} .

The response of the system and its increment can be written in the following matrix form:

17

y_{0} = S A, Δ y = S Δ A .

After substituting Eq. (17) into the incremental Eq. (13) and the unbalanced force Eq. (14), the Galerkin averaging process is applied to obtain the equations for the unknown quantities $Δ A$ and $Δ Ω$ :

18

{\bar{K}}_{j c} Δ A = \bar{R} - {\bar{R}}_{j c} Δ Ω,

where:

{\bar{K}}_{j c} = \int_{0}^{2 π} S^{T} (Ω_{0}^{2} \hat{M} S'' + Ω_{0} \hat{C} S^{'} + \hat{K} S) d τ_{h},

\bar{R} = \int_{0}^{2 π} S^{T} \hat{F} (τ_{h}) - \int_{0}^{2 π} S^{T} (Ω_{0}^{2} \hat{M} S'' + Ω_{0} \hat{C} S^{'} + \hat{K} S)] d τ_{h} A,

{\bar{R}}_{j c} = \int_{0}^{2 π} S^{T} (2 Ω_{0} \hat{M} S'' + \hat{C} S^{'}) d τ_{h} A .

3.2. Arc length extension method

The arc length parameter equation, corresponding to Eq. (11), can be expressed as:

19

g (p) - s = 0 .

Assuming $g (p) = p^{T} p$ and $p = [A^{T}, Ω]^{T}$ and substituting the increments of $A$ , $Ω$ , and $s$ in Eq. (19), the increment equation can be obtained, as shown below:

20

\frac{\partial g}{\partial A^{T}} {Δ A} + \frac{\partial g}{\partial Ω} Δ Ω - Δ s + g - s = 0 .

Fig. 2 shows a part of the analytical balance path of the arc-length extension method.

Eq. (19) can be rewritten as:

21

g (p) - s = {p^{'}}^{T} {p - p_{c}} = 0,

22

P_{u} = P_{c} + Δ s (\frac{\{P_{c} - P_{c c}\}}{‖P_{c} - P_{c c}‖}) .

The initial values of the upper and lower points $P_{u}$ of the balance path are determined by the values of the previous two points, $P_{c}$ , and $P_{c c}$ .

Fig. 2Partial schematic diagram of the balance path based on the arc length extension method

The complete incremental equation can be obtained by combining Eqs. (18) and (20):

23

[{\bar{J}}_{p}] \{Δ p\} = [\begin{matrix} [{\bar{K}}_{j c}] & [{\bar{R}}_{j c}] \\ {\{\frac{\partial g}{\partial A}\}}^{T} & \frac{\partial g}{\partial Ω} \end{matrix}] \{\begin{matrix} ∆ A \\ ∆ Ω \end{matrix}\} = \{\begin{matrix} \bar{R} \\ - g + ∆ s + s \end{matrix}\} = [\tilde{R}],

where $[{\bar{J}}_{p}]$ is the Jacobian matrix relative to ${p}$ .

The above equation is expressed in an iterative form that can be easily calculated, as shown below:

24

p_{i + 1} = p_{i} + [{\bar{J}}_{p}]^{- 1} (p_{i}) \tilde{R} (p_{i}) .

Eq. (24) represents the Newton-Raphson iterative equation that is obtained after introducing the arc length parameter. The arc length parameter is used to predict the value of the next solution from the current solution and is used as the initial value in the next iteration.

4. Application of the improved method

The transmission diagram of the double-wide helical planetary compound system, with a two-level power branch, for high speed and heavy duty applications, is shown in Fig. 3. The system is composed of a star gear train I (sun gear $z_{s}^{Ι}$ , planet gear $z_{m}$ , and ring gear $z_{r}^{Ι}$ ) that is connected to a planet gear train II (sun gear $z_{s}^{Ι Ι}$ , planet gear $z_{n}$ , ring gear $z_{r}^{Ι Ι}$ , and planet carrier $H$ ). The superscripts I and II correspond to the series of the component. The input power is transmitted to the load $L$ by the sun gear $z_{s}^{Ι}$ . The input speed of the second-stage planetary gear train is reduced according to the deceleration of the first stage planetary gear train, thereby increasing the stability and smoothness of the transmission. The system parameters are listed in Table 1 and Table 2. The calculated nonlinear response characteristics of the system are shown in Figs. 4 to 6, which correspond to the time domain response history, phase diagram, and Poincare mapping of the system, respectively. Owing to the space constraints in this manuscript, we have only mentioned the torsional response of the representative components. However, this does not imply that the translational response of the components can be neglected.

Since the number of unknowns exceeds the number of equations, the expected increment must be specified before performing the actual calculation. The increment specified in this article is equal to $Δ Ω$ . Based on the structural flowchart of the improved incremental harmonic balance method depicted in Fig. 4 and the dimensionless parameters listed in Table 1 and Table 2, the specific iterative process of the improved method is given as follows.

Fig. 3Transmission diagram of a double-wide helical planetary compound system with a two-level power branch

(1) The initial value $y_{0}$ is fixed according to the excitation frequency $Ω = Ω_{0}$ .

(2) The increment $Δ A$ is obtained by substituting the value of $Δ Ω$ in Eq. (18). Replace $A$ with $A + Δ A$ to obtain the updated values of $Δ A$ from Eq. (18). This is used to obtain $Δ y$ from Eq. (17). The modified solution $y$ is then obtained from Eq. (12). The process is repeated until the value of $A$ satisfies $R = 0$ .

(3) A new increment is provided to $Ω_{0}$ such that $Ω = Ω_{0} + Δ Ω$ . The value of $A$ that is obtained in (2) is set as the initial value. The harmonic balance process is repeated and the value of $A$ is updated until it meets the condition $R = 0$ .

(4) The arc length parameter $s$ is introduced to determine the initial value of the next point from the value of $A$ and $Ω$ obtained in Eqs. (2) and (3). This is substituted as the initial value in Eq. (18), and the iteration step mentioned in Eq. (2) is repeated.

Fig. 4Flowchart of the numerical integral calculation system response

The amplitude-frequency characteristic curves of the system operating at the approximate working speed are shown in Fig. 5. The figure compares the results obtained by the incremental harmonic balance method and those obtained by the numerical integration (NI) method. The latter is based on the 4th and 5th step variable step Runge-Kutta method.

The meshing frequency lies in the range of 3.3-3.7 kHz, as shown in Figs. 5(a), (b) and (c). Amplitude jumps are observed in both the sun gear and star gear of the star gear system. The two methods are in this area. The amplitude values obtained from both the methods do not match as well in this region as they do in the others. In addition, the results of the numerical integration method indicate the presence of a resonance peak in other regions; however, the IHB method does not seem to indicate the presence of a resonance peak. This is because the number of harmonic response terms is insufficient.

Table 1System parameters

Physical quantity \Transmission parts marking	Sun gear		Ring gear		Planet carrier		Star / Planet gears
Physical quantity \Transmission parts marking	Level 1	Level 2	Level 1	Level 2	Level 1	Level 2	Level 1	Level 2
Quality (kg)	102	398	258	647		1300	300	647
Equivalent moment of inertia $Ι / r^{2}$ (kg)	51.5	199	244	613		1300	150	324
Base circle diameter (mm)	385.27	479.24	1700.84	1700.84			657.78	610.80
Number of teeth	41	51	181	181			70	65
Meshing stiffness (N/m)	$k_{s p}^{Ι} = 21.532 \times 1 0^{9}$ ; $k_{r p}^{Ι} = 3.393 \times 1 0^{9}$ ; $k_{s p}^{Ι Ι} = 21.802 \times 1 0^{9}$ ; $k_{r p}^{Ι Ι} = 3.281 \times 1 0^{9}$
Support stiffness (N/m)	$k_{r}^{Ι} = k_{c}^{Ι Ι} = k_{r}^{Ι Ι} = k_{s}^{Ι Ι} = k_{p}^{Ι} = k_{p}^{Ι Ι} = 1.0 \times 1 0^{10}$ ; $k_{s}^{Ι} = 100$ ; $k_{r s}^{12} = 0$
Torsional stiffness of central member (N/m)	$k_{s u}^{Ι} = k_{r u}^{Ι} = k_{s u}^{Ι Ι} = k_{r u}^{Ι Ι} = k_{c u}^{Ι Ι} = k_{p m u} = k_{p n u} = 1 0^{10}$
Torsional stiffness of shaft (Nm/rad)	$k_{r s u}^{12} = 2.0 \times 1 0^{8}$
Pressure angle ( $°$ )	$α_{s m} = α_{r m} = α_{s n} = α_{r n} = 20$
Helix angle ( $°$ )	$β_{b m} = β_{b n} = 24.62$

Fig. 5Variation of the torsional response amplitude of each component of the planetary transmission system by using the IHB method

a) Ring gear of star gear train

b) Sun gear of star gear train

c) Star gear of star gear train

d) Planet carrier of planetary gear train

e) Sun gear of planetary gear train

f) Planet gear of planetary gear train

Table 2Dimensionless parameters

Parameter representation	Star gear train	Planetary gear train
Dimensionless parameter $\hat{ω}$	$\hat{ω} = 1.1893 e + 004$
Dimensionless parameter ${\bar{e}}_{s 1}^{Ι}$	${\bar{e}}_{s 1}^{Ι} = 3.6059 e - 005$
Dimensionless meshing frequency $ω_{m i}$	$ω_{m 1} = 0.275 8$	$ω_{m 1} = 0.275 8$
Damping ratio $ζ$	$ζ = 0.0 2$
Dimensionless error amplitude	${\hat{e}}_{s m}^{Ι} = {\hat{e}}_{s n}^{Ι Ι} = 1.109 3$

The ring gear of the star gear system, shown in Figs. 5(a), (b) and (c), and the sun gear of the planetary gear system, shown in Fig. 5(e), produces minimal amplitude resonance. The results obtained by the numerical method and IHB method follow a similar trend; however, a significant difference exists between the amplitude values. The variations in the amplitude of the planetary carrier, shown in Fig. 5(d), and the planetary gears, shown in Fig. 5(f), are relatively gentle. The curves obtained from the two methods gradually tend to become consistent with increasing mesh frequency, resulting in a significant amplitude variation of the planetary gears. This is consistent at a 4.0 kHz bit frequency. The above analysis demonstrates that the improved method is in agreement with the law of amplitude-frequency change for each part of the system, thereby illustrating the feasibility of this method.

5. Application of the multi-scale perturbation analysis method

The incremental harmonic balance method, which is based on the arc-length continuation technique, is a semi-analytical and semi-numerical method. It is necessary to perform a purely analytical study of the dynamic characteristic equation of the system. Current studies adopt the multi-scale perturbation analysis method to obtain the analytical solutions of parametric excitation and gap nonlinear system equations [31]. The multi-scale perturbation method (MsPM) can obtain the analytical frequency response functions of a system, including the fundamental, subharmonic, and superharmonic resonance responses. Thus, this technique demonstrates the impact of important parameters on the response of the nonlinear dynamic characteristics, unlike conventional numerical methods.

The small parameter $ε = |{\hat{c}}_{s n}^{(1)}| / k_{s p}^{Ι Ι}$ , ${\hat{c}}_{s n}^{(1)}$ is introduced in the first-order Fourier coefficient of the meshing stiffness of the sun and planetary gears in the planetary gear system. $k_{s p}^{I I}$ refers to the average meshing stiffness and is written in its dimensionless form for each gear pair:

25

k_{s m}^{Ι} (τ) = k_{s p}^{Ι} [1 + ε \sum_{l = 1}^{\infty} (b_{s m}^{(l)} e^{j l ω_{m 1} τ} + c c)],

k_{r m}^{Ι} (τ) = k_{r p}^{Ι} [1 + ε \sum_{l = 1}^{\infty} (b_{r m}^{(l)} e^{j l ω_{m 1} τ} + c c)],

26

k_{s n}^{Ι Ι} (τ) = k_{s p}^{Ι Ι} [1 + ε \sum_{l = 1}^{\infty} (c_{s n}^{(l)} e^{j l ω_{m 2} τ} + c c)],

k_{r n}^{Ι Ι} (τ) = k_{r p}^{Ι Ι} [1 + ε \sum_{l = 1}^{\infty} (c_{r n}^{(l)} e^{j l ω_{m 2} τ} + c c)],

where:

c_{s n}^{(l)} = \frac{{\hat{c}}_{s n}^{(l)}}{|{\hat{c}}_{s n}^{(1)}|} = O (1), c_{r n}^{(l)} = \frac{{\hat{c}}_{r n}^{(l)}}{|{\hat{c}}_{s n}^{(1)}|} \frac{k_{s p}^{Ι Ι}}{k_{r p}^{Ι Ι}} = O (1),

b_{s m}^{(l)} = \frac{k_{s p}^{Ι Ι}}{k_{s p}^{Ι}} \frac{{\hat{b}}_{s m}^{(l)}}{|{\hat{c}}_{s n}^{(1)}|} = O (1), b_{r m}^{(l)} = \frac{k_{s p}^{Ι Ι}}{k_{r p}^{Ι}} \frac{{\hat{b}}_{r m}^{(l)}}{|{\hat{c}}_{s n}^{(1)}|} = O (1) .

The time required for the contact gear pair to disengage is assumed to be negligible with respect to the response period, i.e., $ξ / T = O (ε)$ , where $ξ$ is the disengagement time, and $T = 2 π / Ω$ is the response period. The Fourier expansion of the non-meshed function of the contact gear pair is expressed in terms of the fundamental frequency $Ω$ , as follows:

27

Θ (δ_{s m}^{Ι}) = 1 + ε \sum_{h = 0}^{\infty} ({\hat{g}}_{s m}^{(h)} e^{j h Ω τ} + c c), Θ (δ_{r m}^{Ι}) = 1 + ε \sum_{h = 0}^{\infty} ({\hat{g}}_{r m}^{(h)} e^{j h Ω τ} + c c),

28

Θ (δ_{s n}^{Ι Ι}) = 1 + ε \sum_{h = 0}^{\infty} ({\hat{θ}}_{s n}^{(h)} e^{j h Ω τ} + c c), Θ (δ_{r n}^{Ι Ι}) = 1 + ε \sum_{h = 0}^{\infty} ({\hat{θ}}_{r n}^{(h)} e^{j h Ω τ} + c c) .

The corresponding eigenvalue of Eq. (11) is expressed as:

29

({\hat{K}}_{b} + {\hat{K}}_{0}) V_{i} = c_{i}^{2} \hat{M} V_{i},

where ${\hat{K}}_{0}$ is the linear time-invariant average meshing stiffness matrix, and $\hat{K} (y, t) = {\hat{K}}_{c} + {\hat{K}}_{b} + {\hat{K}}_{0} + {\hat{K}}_{d} (y, t)$ . ${\hat{K}}_{d}$ represents the change range matrix of the average meshing stiffness ${\hat{K}}_{0}$ , with its mean value equal to zero. ${\hat{K}}_{b}$ is the support torsional stiffness matrix (including the support stiffness of the star and planet gears). ${\hat{K}}_{c}$ is the additional stiffness matrix generated due the transient phase angle of the planetary gear. ${\hat{K}}_{c} = {\hat{C}}_{p n}^{Ι Ι} + {\hat{C}}_{s n}^{Ι Ι} + {\hat{C}}_{r n}^{Ι Ι}$ , where ${\hat{C}}_{p n}^{Ι Ι}$ , ${\hat{C}}_{s n}^{Ι Ι}$ , and ${\hat{C}}_{r n}^{Ι Ι}$ are the coefficient matrices related to $C_{p n}^{Ι Ι}$ , $C_{s n}^{Ι Ι}$ , and $C_{r n}^{Ι Ι}$ , respectively.

The vibration mode is given by $V = [V_{1}, \dots, V_{3 (M + N + 4)}]$ , and it satisfies the relation $V^{T} \hat{M} V = I$ . The average stiffness matrix $\hat{K}$ is given below:

30

\hat{K} = {\hat{K}}_{0} + {\hat{K}}_{b} = \sum_{m = 1}^{M} (k_{s p}^{Ι} {\hat{K}}_{s m}^{Ι} + k_{r p}^{Ι} {\hat{K}}_{r m}^{Ι}) + \sum_{n = 1}^{N} (k_{s p}^{Ι Ι} {\hat{K}}_{s n}^{Ι Ι} + k_{r p}^{Ι Ι} {\hat{K}}_{r n}^{Ι Ι}),

where ${\hat{K}}_{s m}^{Ι}$ , ${\hat{K}}_{r m}^{Ι}$ , ${\hat{K}}_{s n}^{Ι Ι}$ , and ${\hat{K}}_{r n}^{Ι Ι}$ are the coefficient matrices related to $k_{s p}^{Ι}$ , $k_{r p}^{Ι}$ , $k_{s p}^{Ι Ι}$ , and $k_{r p}^{I I}$ .

Substituting Eqs. (25)-(28) into Eqs. (11), we obtain:

31

\hat{M} \overset{⃛}{y} + \hat{C} \ddot{y} + {\hat{K}}_{c} y + ({\hat{K}}_{b} + {\hat{K}}_{0}) y + ε \sum_{n = 1}^{N} k_{s p}^{Ι Ι} Q_{s}^{Ι Ι} {\hat{K}}_{s n}^{Ι Ι} y + ε \sum_{n = 1}^{N} k_{r p}^{Ι Ι} Q_{r}^{Ι Ι} {\hat{K}}_{r n}^{Ι Ι} y + ε \sum_{m = 1}^{M} k_{s p}^{Ι} Q_{s}^{Ι} {\hat{K}}_{s m}^{Ι} y + ε \sum_{m = 1}^{M} k_{r p}^{Ι} Q_{r}^{Ι} {\hat{K}}_{r m}^{Ι} y + O (ε^{2}) = \hat{F},

where:

Q_{s}^{Ι} = [\sum_{l = 1}^{\infty} b_{s m}^{(l)} e^{j l ω_{m 1} τ} + \sum_{h = 0}^{\infty} {\hat{g}}_{s m}^{(h)} e^{j h Ω τ}] + c c, Q_{r}^{Ι} = [\sum_{l = 1}^{\infty} b_{r m}^{(l)} e^{j l ω_{m 1} τ} + \sum_{h = 0}^{\infty} {\hat{g}}_{r m}^{(h)} e^{j h Ω τ}] + c c,

Q_{s}^{Ι Ι} = [\sum_{l = 1}^{\infty} c_{s n}^{(l)} e^{j l ω_{m 2} τ} + \sum_{h = 0}^{\infty} {\hat{θ}}_{s n}^{(h)} e^{j h Ω τ}] + c c, Q_{r}^{Ι Ι} = [\sum_{l = 1}^{\infty} c_{r n}^{(l)} e^{j l ω_{m 2} τ} + \sum_{h = 0}^{\infty} {\hat{θ}}_{r n}^{(h)} e^{j h Ω τ}] + c c .

The modal coordinate transformation $y = V z$ can be used to write the additional stiffness matrix coefficients in terms of the small parameters. The resulting modal coordinate form, obtained after the transformation of the Eq. (31), is given as:

32

{\ddot{z}}_{q} + ε λ_{q} {\dot{z}}_{q} + c_{q}^{2} z_{q} + ε n_{c}^{Ι Ι} \sum_{w = 1}^{3 (M + N + 4)} [\sum_{n = 1}^{N} C_{p n}^{Ι Ι} E_{p n q w} + \sum_{n = 1}^{N} C_{s n}^{Ι Ι} E_{s n q w} + \sum_{n = 1}^{N} C_{r n}^{Ι Ι} E_{r n q w}] z_{w}

+ ε \sum_{w = 1}^{3 (M + N + 4)} [\sum_{m = 1}^{M} (k_{s p}^{Ι} Q_{s}^{Ι} G_{s m q w} + k_{r p}^{Ι} Q_{r}^{Ι} G_{r m q w})

+ \sum_{n = 1}^{N} (k_{s p}^{Ι Ι} Q_{s}^{Ι Ι} G_{s n q w} + k_{r p}^{Ι Ι} Q_{r}^{Ι Ι} G_{r n q w})] z_{w} = V_{q}^{T} \hat{F} = f_{q},

where $G_{s m} = V^{T} {\hat{K}}_{s m}^{Ι} V$ , $G_{r m} = V^{T} {\hat{K}}_{r m}^{Ι} V$ , $G_{s n} = V^{T} {\hat{K}}_{s n}^{I I} V$ , $G_{r n} = V^{T} {\hat{K}}_{r n}^{I I} V$ , $E_{p n} = V^{T} {\hat{C}}_{p n}^{I I} V$ , $E_{s n} = V^{T} {\hat{C}}_{s n}^{I I} V$ , $E_{r n} = V^{T} {\hat{C}}_{r n}^{I I} V$ .

$E_{p n q w}$ , $E_{s n q w}$ , and $E_{r n q w}$ represent the elements in the $q$ th row and $w$ th column of the matrices $E_{p n}$ , $E_{s n}$ , and $E_{r n}$ , respectively. $G_{s m q w}$ , $G_{r m q w}$ , $G_{s n q w}$ , and $G_{r n q w}$ represent the elements in the $q^{t h}$ row and the $w$ th column of the matrices $G_{s m}$ , $G_{r m}$ , $G_{s n}$ , and $G_{r n}$ , respectively. The modal damping factor $2 ζ_{q} c_{q}$ has been introduced and rewritten in terms of the small parameter $ε$ as $λ_{q} = 2 ζ_{q} c_{q}$ . $ε n_{c}^{Ι Ι} = ω_{c}^{Ι Ι}$ , where $ω_{c}^{I I}$ is the speed of the second stage planet carrier. The small parameter $ε$ is related to the time change of the phase angle of the planet gear in the second stage.

A multi-scale method is applied by introducing multi-scale variables such as $τ_{n} = ε^{n} τ$ and $z_{q} (τ_{0}, τ_{1}, \dots) = z_{q 0} (τ_{0}, τ_{1}, \dots) + ε z_{q 1} (τ_{0}, τ_{1}, \dots) + O (ε^{2})$ . Based on the abovementioned variables, the perturbation equation with the first approximate solution is proposed, as shown below:

33

\frac{\partial^{2}}{\partial τ_{0}^{2}} z_{q 0} + c_{q}^{2} z_{q 0} = f_{q},

34

\frac{\partial^{2}}{\partial τ_{0}^{2}} z_{q 1} + c_{q}^{2} z_{q 1} = - 2 \frac{\partial^{2} z_{q 0}}{\partial τ_{0} \partial τ_{1}} - λ_{q} \frac{\partial z_{q 0}}{\partial τ_{0}}

- n_{c}^{Ι Ι} \sum_{w = 1}^{3 (M + N + 4)} (\sum_{n = 1}^{N} C_{p n}^{Ι Ι} E_{p n q w} + \sum_{n = 1}^{N} C_{s n}^{Ι Ι} E_{s n q w} + \sum_{n = 1}^{N} C_{r n}^{Ι Ι} E_{r n q w}) z_{w 0}

- \sum_{w = 1}^{3 (M + N + 4)} [\sum_{m = 1}^{M} (k_{s p}^{Ι} Q_{s}^{Ι} G_{s m q w} + k_{r p}^{Ι} Q_{r}^{Ι} G_{r m q w}) + \sum_{n = 1}^{N} (k_{s p}^{Ι Ι} Q_{s}^{Ι Ι} G_{s n q w} + k_{r p}^{Ι Ι} Q_{r}^{Ι Ι} G_{r n q w})] z_{w 0} .

Eq. (34) is the perturbation equation used to calculate the closed solution. The frequency response characteristics of the system under different excitations can be studied by using this equation.

Fig. 6Amplitude-frequency variation curve of each component of the planetary gear system

a) Ring gear of star gear train

b) Sun gear of star gear train

c) Star gear of star gear train

d) Planet carrier of planetary gear train

e) Sun gear of planetary gear train

f) Planet gear of planetary gear train

6. Analysis of the frequency response characteristics

The amplitude-frequency characteristics of the system are obtained and studied after solving the frequency response equations under different resonance conditions according to the multi-scale perturbation analysis method.

A natural frequency of 973.1 Hz is selected to study the frequency response characteristics of the system during resonance in the vibration mode of the planetary gear system. The amplitude-frequency characteristics of the system are analyzed and compared with the results obtained by using the numerical integration method. The calculated amplitude-frequency characteristic curve is shown in Fig. 6. Fig. 6(a) shows that the response amplitude of the ring gear of the star gear train calculated by the multi-scale method is differs significantly from that of the numerical method. However, the trends followed by the variation of the responses in both methods are identical, as shown in Fig. 6(e). The variation of the sun gear in the middle planetary gear system is similar to that of the ring gear in the star gear system. The response amplitudes of the sun gear and the star gear in the star gear system are not consistent with those obtained from the numerical method. The variation trends obtained through both methods were also different. No amplitude jump was observed, as shown in Figs. 6(b) and 6(c). The variation trends of the planetary carrier and gears of the planetary gear train, as shown in Figs. 6(d) and 6(f), are identical. The larger difference is also the magnitude of the amplitude. The difference between the values of the amplitude obtained through both methods is relatively large for the planet carrier and small for the planet gear.

The impact of the variation of the damping ratio on the amplitude-frequency response characteristics, upon the introduction of the 3rd harmonic error, is studied. It can be seen from Fig. 7 that the bifurcation characteristics of the system are complex, and the introduction of errors increases the influence of the damping ratio. After increasing the damping ratio, the response amplitude of sun and star gears in the star gear system is reduced. The component response amplitude has increased. In the follow-up research, the influence of the damping ratio on the system’s amplitude-frequency characteristics can be further confirmed.

Fig. 7Bifurcation diagram of the variation of damping ratio with the error

a) Ring gear of star gear train

b) Sun gear of star gear train

c) Star gear of star gear train

d) Planet carrier of planetary gear train

e) Sun gear of planetary gear train

f) Planet gear of planetary gear train

The steady-state response of the ring and sun gears of the star gear train is either a non-harmonic periodic response or a simple harmonic periodic response, as shown in Figs. 7(a) and (e). The gears of the planetary gear train always maintain a harmonic response without bifurcation, as shown in Fig. 7(f). The steady-state response of the planetary carrier of the planetary gear train, as shown in Fig. 7(d), is bifurcated from a single period to a double period when the damping ratio $ζ$ is equal to 0.03. The steady-state response of the sun gear of the star gear system, as shown in Fig. 7(b), directly branches from period doubling to a harmonic period response at $ζ = 0.035$ , and attains a chaotic state at $ζ = 0.023$ . The steady-state response of the star gear, as shown in Fig. 6(c), involves a period-doubling bifurcation from a single-cycle bifurcation at $ζ = 0.036$ . The steady-state response then bifurcates from period-doubling to a chaotic state at $ζ = 0.023$ .

The nonlinear response characteristics are depicted through a phase diagram of the system-time domain in Fig. 8. The phase diagrams shown in Figs. 8(a), (b), and (c) form a closed curve loop, irregular shape, and an open curve, respectively. The phase diagrams shown in Figs. 8(d) and (f) are ellipses. The response carrier of the planetary gear train produces a pseudo-periodic response, as shown in Fig. 8(e).

The above analysis shows that the results obtained by the multi-scale method can predict the trend of variation of the responses of each component in a few regions. However, it produces linear changes in the region involving an increase in the amplitude. This behavior is attributed to the fact that only one approximate solution was obtained in this study. It is very difficult to increase the order of the solution for a nonlinear system having multiple degrees of freedom. Thus, the numerical method of calculation is more suited to be the main method, with the analytical method serving as an auxiliary option.

Fig. 8Phase plan of the composite transmission system

a) Ring gear of star gear train

b) Sun gear of star gear train

c) Star gear of star gear train

d) Planet carrier of planetary gear train

e) Sun gear of planetary gear train

f) Planet gear of planetary gear train

7. Conclusions

The current paper proposes the application of the semi-numerical incremental harmonic balance and semi-analytical multi-scale perturbation methods to a two-stage series compound planetary transmission system. Through this study, we attempt to solve the dynamic characteristic equation of a two-stage series composite planetary transmission system through an analytical calculation.

The arc-length continuation technology is introduced to improve the incremental harmonic balance method. The improved method is used to calculate and analyze the amplitude-frequency characteristics of the system. The feasibility and effectiveness of the method are verified by comparing the results with those obtained using the numerical integration method.

The analytical multi-scale perturbation method is then applied to a two-stage series compound planetary transmission system. The frequencies of the main resonance, subharmonic resonance, and superharmonic resonance are obtained. A comparison between the current results and those obtained from the numerical integration method suggests that it is feasible to use the multi-scale method to analyze the two-stage series compound planetary transmission system. However, the results may not be accurate in some regions.

References

Inalpolat M., Kahraman A. A theoretical and experimental investigation of modulation sidebands of planetary gear sets. Journal of Sound and Vibration, Vol. 323, Issues 3-5, 2009, p. 677-696.

Publisher
Feng Z. P., Zuo M. J. Vibration signal models for fault diagnosis of planetary gearboxes. Journal of Sound and Vibration, Vol. 331, Issue 22, 2012, p. 4919-4939.

Publisher
Hotait M. A., Kahraman A. Experiments on the relationship between the dynamic transmission error and the dynamic stress factor of spur gear pairs. Mechanism and Machine Theory, Vol. 70, 2013, p. 116-128.

Publisher
He G. L., Ding K., Li W. H., Li Y. Z. Frequency response model and mechanism for wind turbine planetary gear train vibration analysis. IET Renewable Power Generation, Vol. 11, Issue 4, 2017, p. 425-432.

Publisher
Pan C. D., Yu L., Liu H. L., Chen Z. P., Luo W. F. Moving force identification based on redundant concatenated dictionary and weighted l₁-norm regularization. Mechanical Systems and Signal Processing, Vol. 98, 2018, p. 32-49.

Publisher
Chen Z., Chan T. H. T., Nguyen A., Yu L. Identification of vehicle axle loads from bridge responses using preconditioned least square QR-factorization algorithm. Mechanical Systems and Signal Processing, Vol. 128, 2019, p. 479-496.

Publisher
Suslin A., Pilla C. Study of loading in point-involute gears. Procedia Engineering, Vol. 176, 2017, p. 12-18.

Publisher
Morgado T. L. M., Branco C. M., Infante V. A failure study of housing of the gearboxes of series 2600 locomotives of the Portuguese Railway Company. Engineering Failure Analysis, Vol. 15, Issues 1-2, 2008, p. 154-164.

Publisher
Hu W. G., Liu Z. M., Liu D. K., Hai X. Fatigue failure analysis of high speed train gearbox housings. Engineering Failure Analysis, Vol. 73, 2017, p. 57-71.

Publisher
Weis P., Kučera Ľ., Pecháč P., Močilan M. Modal Analysis of Gearbox Housing with Applied Load. Procedia Engineering, Vol. 192, 2017, p. 953-958.

Publisher
Lin T. L., Zhang K. An analytical study of the free and forced vibration response of a ribbed plate with free boundary conditions. Journal of Sound and Vibration, Vol. 422, 2018, p. 15-33.

Publisher
Zhou H. A., Zhao Y. G., Wu H. Y., Meng J. B. The vibroacoustic analysis of periodic structure-stiffened plates. Journal of Sound and Vibration, Vol. 481, 2020, p. 115402.

Publisher
Sánchez M. B., Pleguezuelos M., Pedrero J. I. Approximate equations for the meshing stiffness and the load sharing ratio of spur gears including hertzian effects. Mechanism and Machine Theory, Vol. 109, 2017, p. 231-249.

Publisher
Liang X., Zuo M. J., Feng Z. Dynamic modeling of gearbox faults: a review. Mechanical Systems and Signal Processing, Vol. 98, 2018, p. 852-876.

Publisher
Ege K., Roozen N. B., Leclère Q., Rinaldi R. G. Assessment of the apparent bending stiffness and damping of multilayer plates; modelling and experiment. Journal of Sound and Vibration, Vol. 426, 2018, p. 129-149.

Publisher
Marchetti F., Ege K., Leclère Q., Roozen N. B. On the structural dynamics of laminated composite plates and sandwich structures; a new perspective on damping identification. Journal of Sound and Vibration, Vol. 474, 2020, p. 115256.

Publisher
Garambois P., Donnard G., Rigaud E., Perret Liaudet J. Multiphysics coupling between periodic gear mesh excitation and input/output fluctuating torques: Application to a roots vacuum pump. Journal of Sound and Vibration, Vol. 405, 2019, p. 158-174.

Publisher
Acri A., Nijman E., Conrado E., Offner G. Experimental structure-borne energy flow contribution analysis for vibro-acoustic source ranking. Mechanical Systems and Signal Processing, Vol. 115, 2019, p. 753-768.

Publisher
Yang Y., Fenemore C., Kingan M. J., Mace B. R. Analysis of the vibroacoustic characteristics of cross laminated timber panels using a wave and finite element method. Journal of Sound and Vibration, Vol. 494, 2021, p. 115842.

Publisher
Wu H., Wu P. B., Li F. S., Shi H. L., Xu K. Fatigue analysis of the gearbox housing in high-speed trains under wheel polygonization using a multibody dynamics algorithm. Engineering Failure Analysis, Vol. 100, 2019, p. 351-364.

Publisher
Wang Q. B., Zhao B., Fu Y., Kong X. G., Ma H. An improved time-varying mesh stiffness model for helical gear pairs considering axial mesh force component. Mechanical Systems and Signal Processing, Vol. 106, 2018, p. 413-429.

Publisher
Dai H., Long X. H., Chen F., Bian J. Experimental investigation of the ring-planet gear meshing forces identification. Journal of Sound and Vibration, Vol. 493, 2021, p. 115844.

Publisher
Kosała K. Calculation models for analysing the sound insulating properties of homogeneous single baffles used in vibroacoustic protection. Applied Acoustics, Vol. 146, 2019, p. 108-117.

Publisher
Rosa S. D., Desmet W., Ichchou M., Ouisse, Scarpa F. Vibroacoustics of periodic media: multi-scale modelling and design of structures with improved vibroacoustic performance. Mechanical Systems and Signal Processing, Vol. 142, 2020, p. 106870.

Search CrossRef
Bi S. F., Ouisse M., Foltête E., Jund A. Virtual decoupling of vibroacoustical systems. Journal of Sound and Vibration, Vol. 401, 2017, p. 169-189.

Publisher
Arasan U., Marchetti F., Chevillotte F., Tanner G., Chronopoulos D., Gourdon E. On the accuracy limits of plate theories for vibro-acoustic predictions. Journal of Sound and Vibration, Vol. 493, 2021, p. 115848.

Publisher
Daneshjou K., Talebitooti R., Kornokar M. Vibroacoustic study on a multilayered functionally graded cylindrical shell with poroelastic core and bonded-unbonded configuration. Journal of Sound and Vibration, Vol. 393, 2017, p. 157-175.

Publisher
Rohan E., Lukeš V. Homogenization of the vibro-acoustic transmission on perforated plates. Applied Mathematics and Computation, Vol. 361, 2019, p. 821-845.

Publisher
Guo Y., Eritenel T., Ericson T. M., Parker R. G. Vibro-acoustic propagation of gear dynamics in a gear-bearing-housing system. Journal of Sound and Vibration, Vol. 333, Issue 22, 2014, p. 5762-5785.

Publisher
Tomilina T. M. New approaches to design of structures with required vibroacoustic properties. Procedia Engineering, Vol. 106, 2015, p. 350-353.

Publisher
Li Y. Z., Ding K., He G. L., Yang X. Q. Vibration modulation sidebands mechanisms of equally-spaced planetary gear train with a floating sun gear. Mechanical Systems and Signal Processing, Vol. 129, 2019, p. 70-90.

Publisher

About this article

Received

06 December 2020

Accepted

21 March 2021

Published

06 May 2021

SUBJECTS

Chaos, nonlinear dynamics and applications

DOI

https://doi.org/10.21595/jve.2021.21832

Keywords

two-stage series composite system

dual-power branches

planetary transmission

improved incremental harmonic balance method

multi-scale perturbation method

dynamic characteristics analysis

Acknowledgements

The authors would like to thank the Northeast Forestry University (NEFU), Heilongjiang Institute of Technology (HLJIT), and the Harbin Institute of Technology (HIT) for their support.

The research topic was supported by the Doctoral Research Startup Foundation Project of Heilongjiang Institute of Technology (Grant No. 2020BJ06, Yongmei Wang, HLJIT), the Natural Science Foundation Project of Heilongjiang Province (Grant No. LH2019E114, Baixue Fu, HLJIT), the Basic Scientific Research Business Expenses (Innovation Team Category) Project of Heilongjiang Institute of Engineering (Grant No. 2020CX02, Baixue Fu, HLJIT), and the Special Project for Double First-Class-Cultivation of Innovative Talents (Grant No. 000/41113102, Jiafu Ruan, NEFU).

This is an open access article distributed under the Creative Commons Attribution License, which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.

Previous article in issue Previous Next article in issue Next

Research article

2023 10 31

Dynamic characteristics analysis for gear transmission system in shearer cutting section under different loads

Guo Ye, Xing Deng, Jinyong Ju, Lianchao Sheng

Research article

2019 09 30

Nonlinear dynamics of planetary gear wear in multistage gear transmission system

Xin Wang

Research article

2017 09 30

Analytical coupling characterization of multi-stage planetary gear free vibration considering flexible structure

Aiqiang Zhang, Jing Wei, Datong Qin, Shaoshuai Hou

Research article

2016 11 15

Nonlinear dynamic analysis of a coupled lateral-torsional spur gear with eccentricity

Zhaohui Ren, Jialin Li, Kai Wang, Shihua Zhou

Y. Wang, X. Wang, B. Fu, J. Ruan, and S. An, “Semi-numerical analysis of a two-stage series composite planetary transmission considering IHB and MsP methods,” Journal of Vibroengineering, Vol. 23, No. 3, pp. 708–727, May 2021, https://doi.org/10.21595/jve.2021.21832

Copy Extrica

Copied to clipboard!

TY  - JOUR
DO  - 10.21595/jve.2021.21832
UR  - https://doi.org/10.21595/jve.2021.21832
TI  - Semi-numerical analysis of a two-stage series composite planetary transmission considering IHB and MsP methods
T2  - Journal of Vibroengineering
AU  - Wang, Xigui
AU  - Wang, Yongmei
AU  - Fu, Baixue
AU  - Ruan, Jiafu
AU  - An, Siyuan
PY  - 2021
DA  - 2021/05/06
PB  - JVE International Ltd.
SP  - 708-727
IS  - 3
VL  - 23
SN  - 1392-8716
SN  - 2538-8460
ER  - 

Copy Ris

Copied to clipboard!

@article{Wang_2021,
	doi = {10.21595/jve.2021.21832},
	url = {https://doi.org/10.21595/jve.2021.21832},
	year = 2021,
	month = {may},
	publisher = {{JVE} International Ltd.},
	volume = {23},
	number = {3},
	pages = {708--727},
	author = {Xigui Wang and Yongmei Wang and Baixue Fu and Jiafu Ruan and Siyuan An},
	title = {Semi-numerical analysis of a two-stage series composite planetary transmission considering {IHB} and {MsP} methods},
	journal = {Journal of Vibroengineering}
}

Copy Bibtex

Copied to clipboard!

[1]X. Wang, Y. Wang, B. Fu, J. Ruan, and S. An, “Semi-numerical analysis of a two-stage series composite planetary transmission considering IHB and MsP methods,” Journal of Vibroengineering, vol. 23, no. 3, pp. 708–727, May 2021, doi: 10.21595/jve.2021.21832.

Copy IEEE

Copied to clipboard!

Wang, Xigui, Yongmei Wang, Baixue Fu, Jiafu Ruan, and Siyuan An. “Semi-Numerical Analysis of a Two-Stage Series Composite Planetary Transmission Considering IHB and MsP Methods.” Journal of Vibroengineering 23, no. 3 (May 6, 2021): 708–27. https://doi.org/10.21595/jve.2021.21832.

Copy Chicago

Copied to clipboard!

Semi-numerical analysis of a two-stage series composite planetary transmission considering IHB and MsP methods

Abstract

Highlights

1. Introduction

2. Mathematical model of the modal analysis

2.1. Dimensionless differential equation of the second-stage star gear train

2.2. Dimensionless differential equation of the second-stage planetary gear train

3. Incremental harmonic balance and arc length extension method

3.1. Incremental harmonic balance method

3.2. Arc length extension method

4. Application of the improved method

5. Application of the multi-scale perturbation analysis method

6. Analysis of the frequency response characteristics

7. Conclusions

References

About this article

Related Articles