Semi-analytic semi-numerical analysis of dynamic characteristics of a two-stage coupled series PGT based on IHB and MsP methods

This paper conducts an analytical investigation of the dynamic response characteristics of a two-stage coupled series system (TsCSS) with a planetary gear transmission (PGT), 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, in order to verify the feasibility and efficiency of the improved method. Afterwards, the multi-scale perturbation (MsP) method is applied to the proposed TsCSS. The frequency response equations of the primary resonance, subharmonic resonance and superharmonic resonance, are solved in order to generate the frequency response characteristic curves of the planetary transmission system. A comparison between the results obtained by the MsP method and the numerical integration method demonstrates that the former is ideal and credible in most of the regions. Considering the parameters of TsCSS excitation frequency and damping, the nonlinear response characteristics of steady-state motion are mutually converted. The impact of the time-varying parameters and nonlinear deenthing caused by the gear teeth clearance on the amplitude–frequency characteristics of the TsCSS components, is finally studied.

applications of several planetary gear mechanisms and their transmission systems [3][4][5]. Several researchers tackled the simulation analysis methods, in order to gain an in-depth understanding of the dynamic meshing transients of planetary gears [6][7][8]. Several studies were conducted on the dynamic behavior analysis method, which plays an important role in determining the performance of planetary gears. In addition, the machinery industry focuses on the dynamic characteristics of the PGT and the resulting vibrations and noise [9][10][11].
A two-stage series composite system (TsCSS), which consists of a PGT system with dual-power branches, is a complex and flexible mechanical system comprising several parts [12]. A TsCSS is mainly 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]. Considering the characteristics of the type and structure of various acoustic excitation incident fields, some descriptions of the solution process (analytical, experimental and numerical) are proposed, the review is centralized through other matters such as control and optimization techniques in order to improve the vibroacoustic behavior of structures [14]. An analytical strategy is proposed through the double Fourier series solution to acoustically determine the vibration response of a doubly curved composite shell in a thermal environment, revealing the effect of thermal loads on the sound transmission loss (STL). Based on Hamilton's principle and considering shear deformation shallow shell theory (SDSST), the motion equations of the structure integrated with thermal stresses are derived [15]. The increasing number of studies tackling the vibrations and noise produced by underwater devices, particularly on controlling the noise produced by the power rear transmission system of underwater devices, should account for the TsCSS and its meshing process, in addition to the real-time dynamic meshing force acting on the system. Few theoretical studies consider the TsCSS with a PGT, which consists of dual-power branches [16]. It is crucial to possess some professional design knowledge, when studying the unique kinematics and geometric characteristics of a TsCSS consisting of a PGT with dual-power branches [17]. The TsCSS with a dual-power branch PGT, is preferred to the horizontal shaft gear deceleration system, especially in the applications that require high linear speed power density designs and kinematic flexibility, in order to optimize different speed ratios [18][19][20]. It has been demonstrated that reducing the spoke thickness to increase gear flexibility, also resolves several internal gear, planetary frame and operational errors. In addition, it makes the system lighter [21]. Moreover, a flexible internal gear improves the load sharing between planets, which is a crucial feature if the manufacturing and assembly related gear and carrier errors are inevitable. Thus, it is difficult to quantify the factors that influence the TsCSS with a dual-power branch PGT, under quasi-static conditions.
Modals are the inherent characteristics of gear transmission systems [22][23][24]. 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. It is also the starting point for the subsequent transient dynamics, harmonic response and spectrum analyses [25]. Each mode has a corresponding natural frequency, damping ratio and mode shape [26]. A modal analysis is a modern technique used to study the dynamic characteristics of a transmission system structure, including the modal analysis of linear vibration theory and experimental modal analysis [27]. An experimental modal analysis can only be performed after the components of the structure have been processed and assembled. However, it has a longer test cycle and it 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 [28]. Consequently, the experimental analysis is used to verify the results of the theoretical model, and modify it accordingly. Several modal parameters, such as the modal frequency, shape, quality, stiffness and damping, affect the dynamic load design. Considering that the wellknown approximate analytical methods are generally difficult to deal with the periodic orbit bifurcation problem of the dynamic system. At present, the methods for studying the bifurcation of periodic orbits and chaos of dynamic systems mainly include numerical methods, experimental methods and their combinations. These research methods have been applied to nonlinear vibrator systems, axial motion, mechanical vibration systems and other nonlinear dynamic systems, and partially revealed a large number of nonlinear dynamic phenomena, including the bifurcation of infinite loop multiplication period orbits, and Neimark-Sacker bifurcation into chaos etc. Throughout the above literature, semi-numerical semi-analytical methods are still relatively rare.
Due to the fact that most of the internal, planetary and sun gears are excluded from these models, it is impossible to study the impact of the internal gear thickness on the performance of the TsCSS and the stress acting on the planetary, sun gears and load sharing between the planets. Similarly, it is impossible to accurately predict the shape and deflection of the gears. Although these studies demonstrate 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 result in increasing the gear contact stress. These static analyses cannot solely predict the actual design of the system, because the increasing flexibility of the TsCSS causes its performance to change only under dynamic conditions. This may also lead to an increase in the stress acting on the gear. Therefore, this paper studies the inherent characteristics of the TsCSS, constructs its dynamic model, and synthetically analyzes its inherent and dynamic response characteristics.

Mathematical model of the semi-numerical analysis
The first step of the dynamic calculation consists in determining 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 composite transmission system test bench is considered as a single elastic body. The differential equation of motion describing the overall vibration of the test bench is given by: whereü(t) ,u(t) and u(t) respectively denote the acceleration, velocity and displacement vectors of the nodes in the vibration system, M, C and K respectively represent the mass, damping and stiffness matrices of the vibration system, and Q(t) is the external force vector received by the 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 expressed as: The resonance solution form of the equation is given by: where u is the displacement vector, φ is the characteristic vector of amplitude of the displacement vector u, and ω n is the natural angular frequency. The resonance form of the equation is the key to the numerical solution. It is derived while assuming that all the degrees of freedom of the vibrating structure move in a synchronous manner. During this process, the basic shape of the structure does not change, while only the amplitude varies. According to the dimensionless differential equation of the planetary gear train, the matrix form is rewritten as: It is important to mention that the stiffness matrix K is no longer symmetric, due to the transient nature of the phase angle of the planetary gears.

Incremental harmonic balance and arc length extension method
The harmonic balance method uses a description function in order to approximate the nonlinearity caused by the gap. It is widely used in PGT systems. The excitation and response parameters are considered as harmonic functions, substituted in the nonlinear equation. Therefore, the approximate expressions of the response and phase parameters can be obtained 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, due 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, which improves the accuracy of the existing harmonic balance method. A Taylor series expansion is performed on the nonlinear differential equations while ignoring the higher-order derivatives, in order to obtain the differential equations in an incremental form. The Fourier series and the Galerkin method are then used to obtain the 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. Therefore, it is efficient for solving complex nonlinear problems [29]. Rohan and Lukeš used the incremental harmonic balance method for a two-stage star gear train with multiple degrees of freedom. However, the response is only used as an incremental parameter [30]. 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, in order to smoothly track the path through the singular point [31]. The harmonic balance method is based on the continuation of the arc-length. It has been applied to the pure torsion model of a fixed shaft, and a single-stage PGT [32]. The incremental harmonic balance method is based on the continuation of the arc length. It is used to calculate the dynamic characteristic equation of the system used in this paper. The two incremental parameters (response and fundamental frequency) are expressed in terms of the arc length, in order to smoothly overcome the singular points along the path. The formula allowing to calculate the steady-state response of a two-stage herringbone PGT system, is presented in the sequel. Note that this method has not yet been applied to the PGT systems.

Incremental harmonic balance method
In this method, a new time variable τ h τ is introduced. The expression of y, initially written in terms of τ , is rewritten in terms of τ h . Consequently, Eq. (4) is written as: The incremental process is the first component of the IHB method. If y j0 and 0 are the solutions of Eq. (5), their neighboring points can be expressed as: where y j and are the incremental parameters. By substituting Eq. (6) into Eq. (5) and omitting the high-order small quantities, the incremental equation matrix can be obtained with y j and as the unknown quantities: where R is the unbalanced force vector, also known as the residual correction term.
If y j0 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 it is given by: The response of the system and its increment can be written in the following matrix form: After substituting Eq. (11) into the incremental Eq. (7) and the unbalanced force Eq. (8), the Galerkin averaging process is performed to obtain the equations of the unknown quantities A and where The arc length parameter equation, corresponding to Eq. (5), can be expressed as: The increment equation can be obtained by assuming that g( p) p T p, p [A T , ] T and substituting the increments of A, and S in Eq. (13): Figure 1 shows a part of the analytical balance path of the arc-length extension method. Equation (13) can be rewritten as: 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 cc .
The complete incremental equation can be obtained by combining Eqs. (12) and (14): where J p is the Jacobian matrix relative to { p}. Equation (17) is expressed in an iterative form that can be easily computed as: Equation (18) represents the Newton-Raphson iterative equation, which 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. In addition, it is used as the initial value in the next iteration.

Analysis application of the improved method
The transmission diagram of the double-wide helical planetary composite system, with a two-level power branch, for high speed and heavy-duty applications, is shown in Fig. 2. The system is composed of a star gear train I (sun gear z I s , planet gear z m and ring gear z I r ), that is connected to a planet gear train II (sun gear z II s , planet gear z n , ring gear z II r and planet carrier H). The superscripts I and II correspond to the series of the    Tables 1 and 2. The calculated nonlinear response characteristics of the system are shown in Figs. 3, 4, 5, that correspond to the time domain response history, phase diagram and Poincare mapping of the system, respectively. Due to the space constraints in this paper, the torsional response of the representative components, is only mentioned. 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 should be specified before performing the calculation. In this paper, an increment of is used. Based on the structural flowchart of the improved incremental harmonic balance method (cf. Figure 3) and the dimensionless parameters (cf. Tables 1 and 2), the specific iterative process of the improved method is summarized as follows.
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). A is replaced with A+ A in order to obtain the updated values of parameter y 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 obtained in Eq. (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 values of A and obtained in Eqs. (2) and (3), respectively. This is substituted as the initial value in Eq. (18), and the iteration step in Eq. (2) is repeated.
The amplitude-frequency characteristic curves of the system, operating at the approximate working speed, are presented in Fig. 4. 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 Fig. 4b and c. Amplitude jumps are observed in the sun gear and star gear, of the star gear system. Therefore, the two methods are in this area. The amplitude values obtained from the two methods do not match in this region, as well as in the other regions. In addition, the results of the numerical integration method demonstrate the presence of a resonance peak in other regions. However, the IHB method does not indicate the presence of a resonance peak. This is due to the fact that the number of harmonic response terms is insufficient.
The ring gear of the star gear system (cf. Fig. 4a), and the sun gear of the planetary gear system (cf. Fig. 4e), produce a minimal amplitude resonance. The results obtained by the numerical method and IHB method have a similar trend. However, a significant difference exists between the amplitude values. The variations in the amplitude of the planetary carrier (cf. Fig. 4d) and the planetary gears (cf. Fig. 4f), are relatively moderate. The curves obtained from the two methods gradually tend to become consistent with the increase of the mesh frequency, which results 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 coherent with the law of amplitude-frequency change for each part of the system, which demonstrates the feasibility of this method.

Simulation application of the multi-scale perturbation analysis method
Based on the time-varying nature of the meshing stiffness of the gear pair and the disengagement characteristics of the gear teeth in the PGT system, the dynamic equation is a multi-degree-of-freedom system of nonlinear differential equations, which is a destined result. The incremental harmonic balance method, which is based on the arc-length continuation technique, is a semi-analytical and semi-numerical method. In fact, it is necessary to perform a purely analytical study of the dynamic characteristic equation of the system. The current studies use the multi-scale perturbation analysis method to obtain the analytical solutions of parametric excitation, and gap nonlinear system equations [33]. The multi-scale perturbation method (MsPM) can obtain the analytical frequency response functions of a system, including the fundamental, subharmonic and superharmonic resonance responses [34]. Thus, this technique demonstrates the impact of the crucial parameters on the response of the nonlinear dynamic characteristics, in contrast to the conventional numerical methods.
The small parameter ε ĉ sn is introduced into the first-order Fourier coefficient of the sun gear and planetary gears meshing stiffness, in the planetary gear system. k II sp refers to the average meshing stiffness. It is written in its dimensionless form for each gear pair as: where c 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 : The corresponding eigenvalue of Eq. (11) is expressed as: whereK 0 is the linear time-invariant average meshing stiffness matrix,K (y, t) K e +K b +K 0 +K d (y, t), K d represents the change range matrix of the average meshing stiffnessK 0 with zero mean,K b is the support torsional stiffness matrix (including the support stiffness of the star and planet gears),K e is the additional stiffness matrix generated due the transient phase angle of the planetary gear,K e Ĉ II pn +Ĉ II sn +Ĉ II rn wherê C II pn ,Ĉ II sn andĈ II rn are the coefficient matrices related to C II pn , C II sn and C II rn , respectively. The vibration mode is given by V [V 1 , · · · , V 3(M+N +4) ]. It satisfies the relation V TM V I . The average stiffness matrix is given by: 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 Eq. (25), is given by: where , E pnqw ,E snqw ,and E rnqw respectively represent the elements in the q th row and ω th column of matrices E pn ,E sn and E rn , G smqw ,G rmqw ,G snqw and G rnqw represent the elements in the q th row and ω th column of matrices G sm ,G rm ,G sn and G rn , 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 . In addition, εn II c ω II c , where ω II c is the speed of the second stage planet carrier. The small parameter ε is related to the time change of the planet gear phase angle in the second stage. A multiscale method is performed by introducing multi-scale variables such as τ n ε n τ and z q (τ 0 , τ 1 , · · ·) z q0 (τ 0 , τ 1 , · · ·) + εz q1 (τ 0 , τ 1 , · · ·) + O(ε 2 ), for example. Based on these variables, the perturbation equation with the first approximate solution is proposed: Equation (28) is the perturbation equation used to calculate the closed solution. Using this equation, the frequency response characteristics of the system can be studied under different excitations.
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 using the numerical integration method. The calculated amplitude-frequency characteristic curve is presented in Fig. 5.
It can be seen from Fig. 5a that the response amplitude of the ring gear of the star gear train, calculated using the multi-scale method, highly differs from that of the numerical method. In addition, no amplitude jump is observed, as shown in Fig. 5b and c. The variation trends of the planetary carrier and gears of the planetary gear train are identical, as shown in Fig. 5d and f. The trends followed by the variation of the responses in the two methods are identical, as shown in Fig. 5e. The magnitude of the amplitude is also the larger difference. The difference between the amplitude values obtained by the two methods is relatively large for the planet carrier, and small for the planet gear.

Validation based on frequency response characteristic analysis
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 by the two methods are also different. 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. 6 that the bifurcation characteristics of the system are complex, and the introduction of errors increases the influence of the damping ratio. In addition, 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 Fig. 6a and e. The gears of the planetary gear train always maintain a harmonic response without bifurcation, as shown in Fig. 6f. Moreover, the steady-state response of the planetary carrier of the planetary gear train (cf. Fig. 6d), 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 (cf. Fig. 6b), directly branches from a 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 (cf. Fig. 6c), 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. Furthermore, the nonlinear response characteristics are illustrated by a phase diagram of the system-time domain in Fig. 7. The phase diagrams shown in Fig. 7a, b, and c form a closed curve loop, irregular shape and open curve, respectively. The phase diagrams shown in Fig. 7d and f are ellipses. The response carrier of the planetary gear train produces a pseudo-periodic response, as shown in Fig. 7e.
This analysis demonstrates that the results obtained by the multi-scale method, can predict the variation trend of the responses of each component in few regions. However, it produces linear changes in the region, which involves an amplitude increase. This behavior is due to the fact that only one approximate solution is 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 calculation method should be considered as the main method, while the analytical method can serve as an auxiliary option.

Conclusion
This paper proposes the application of the semi-numerical incremental harmonic balance and semi-analytical multi-scale perturbation methods, to a two-stage series composite planetary transmission system. Using this study, the dynamic characteristic equation of a two-stage series composite planetary transmission system, is solved by analytical calculations.
The arc-length continuation technology is introduced to improve the incremental harmonic balance method. The improved method is then used to calculate and analyze the amplitude-frequency characteristics of the system. The feasibility and efficiency of the method are verified by comparing its results with those obtained using the numerical integration method.
Afterwards, the analytical multi-scale perturbation method is applied to a two-stage series composite planetary transmission system. The frequencies of the main resonance, subharmonic resonance and superharmonic resonance are obtained. Finally, a comparison between the current results and those obtained from the numerical integration method, demonstrates that the multi-scale method can be used to analyze the two-stage series composite planetary transmission system. However, the results may not be accurate in some regions. This study performs an analytical study of the dynamic response characteristics of a two-stage series composite system with a PGT system that consists of dual-power branches. We believe that our study makes a significant contribution to the literature because it provides an alternative to the computationally heavy and financially expensive numerical methods to analyze PGT systems. It provides a platform for researchers to build upon this study and develop a purely analytical technique with higher accuracy. Further, we believe that this paper will be of interest to the readership of your journal because the frequency response equations of the primary resonance, subharmonic resonance, and superharmonic resonance are solved to generate the frequency response characteristic curves of the PGT 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 TsCSS 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 TsCSS components are studied in this special topic.
Our main contribution is worth mentioning that this paper is a sub project of China Naval Facilities National Defense Construction Project, which is a part research content of the transmission gears thermal deformation topic of warship power rear drive system, these research results have been applied in warship power rear drive system.

Contribution of each individual co-author
No conflict of interest exits in the submission of this manuscript, and manuscript is approved by all authors for publication. We would like to declare on behalf of our co-authors that the work described was original research that has not been published previously, and not under consideration for publication elsewhere, in whole or in part. All the authors listed have approved the manuscript that is enclosed.