Nonlinear Frequency Response Analysis of Double-helical Gear Pair Based on the Incremental Harmonic Balance Method

+e torsional dynamic model of double-helical gear pair considering time-varying meshing stiffness, constant backlash, dynamic backlash, static transmission error, and external dynamic excitation was established.+e frequency response characteristics of the system under constant and dynamic backlashes were solved by the incremental harmonic balance method, and the results were further verified by the numerical integration method. At the same time, the influence of time-varying meshing stiffness, damping, static transmission error, and external load excitation on the amplitude frequency characteristics of the system was analyzed. +e results show that there is not only main harmonic response but also superharmonic response in the system. +e time-varying meshing stiffness and static transmission error can stimulate the amplitude frequency response of the system, while the damping can restrain the amplitude frequency response of the system. Changing the external load excitation has little effect on the amplitude frequency response state change of the system. Compared with the constant backlash, increasing the dynamic backlash amplitude can further control the nonlinear vibration of the gear system.


Introduction
Double-helical gear is widely used in vehicle, ship, and other heavy machinery transmission systems because of its strong bearing capacity, compact structure, and axial force cancellation. Its nonlinear vibration characteristics will directly affect the stability and reliability of the transmission system. Comparin and Singh [1] established a three-degree of freedom spur gear transmission system considering the random disturbance of low frequency external excitation caused by torque fluctuation and studied the influence of different backlash, gear damping ratio, and random disturbance on the dynamic behavior of the gear transmission system. Kahraman and Singh [2][3][4] analyzed the nonlinear dynamic frequency response characteristics of spur gear pair under internal excitation by the numerical method and analytical method and found that there are many nonlinear phenomena such as the jump phenomenon and multivalue solution in the vibration response of spur gear. Shen et al. [5,6] established the torsional vibration model of spur gear pair and gave the unified form of nonlinear analytical solution of the gear system including time-varying meshing stiffness and backlash. Wang et al. [7] established the nonlinear dynamic model of spur gear pair with fault and obtained the phase diagram and bifurcation diagram of the single-stage gear system but did not study the nonlinear amplitude frequency characteristics of the gear system with or without faults. Li et al. [8] established the pure torsion dynamic model of spur gear pair and studied the influence of system parameters on the nonlinear dynamic characteristics of the spur gear system. Zhao et al. [9] analyzed the influence of gear side backlash, eccentricity, and tooth surface friction on the nonlinear vibration of the high-speed gearbox system of wind turbine. Sun et al. [10] developed the application of the single degree of freedom analytical harmonic balance method, derived the calculation formula of the analytic harmonic balance method for differential equations of the planetary gear system, and obtained the nonlinear characteristics of the planetary gear transmission system. Yuan et al. [11] studied the influence of cumulative pitch error on the dynamic characteristics of the double-helical gear system. Wang et al. [12,13] studied the vibration characteristics of the double-helical gear transmission system under different loads and different coincidence degrees. In reference [9,[11][12][13], the vibration problem of the gear system is studied, but the strong nonlinear characteristics caused by backlash are not further explained. In fact, it still belongs to the category of linear dynamics. In addition, many scholars have performed a lot of research on the load sharing characteristics of double-helical gears [14][15][16][17]. Mo et al. [18] analyzed the nonlinear dynamic characteristics of double frequency problem of the double-helical gear closed differential planetary transmission system under no-load or light load conditions. Shin and Palazzolo [19] proposed a new method for modeling and analyzing the geared rotorbearing system. Dong et al. [20] established a 6-DOF spur gear dynamic model including time-varying meshing stiffness, backlash, support backlash, and tooth surface friction. Yavuz et al. [21] first proposed a nonlinear dynamic model of parallel shaft gear and intersecting shaft gear and solved the nonlinear algebraic equations by the multiorder harmonic balance method combined with continuous Fourier transform, but the model is complex and not widely applicable. Xiang et al. [22] established a torsional nonlinear dynamic model of multistage planetary gear considering time-varying meshing stiffness, comprehensive gear error, and piecewise backlash nonlinearities and studied the effects of excitation frequency, backlash, and damping on bifurcation characteristics. Yang et al. [23] established a single degree of freedom spur gear pair model including time-varying mesh stiffness and static transmission error. e nonlinear frequency response characteristics of a spur gear system were studied. Lian et al. [24] established a multiple degrees of freedom nonlinear dynamic model of a gear pair with timevarying mesh stiffness, mesh damping, backlash, dynamic transmission errors, and radial clearance of ball bearing by using the mass centralized method. e effects of main parameters such as exciting frequency, radial clearance of ball bearing, and meshing stiffness ratio on chaos and bifurcation of the system were studied. Li and Peng [25] proposed a gear nonlinear model of multidegree of freedom considering the characteristics of gear tooth. Fractal theory was used to calculate the nonlinear backlash from the tribological aspect. Li et al. [26] studied the influence of the fractal backlash on the dynamic characteristics of the gear bearing system, and the fractal behavior of the backlash was described by W-M function. However, the influence of fractal on the nonlinear frequency response of gears is not studied in reference [25,26].
To sum up, the research on the nonlinear characteristics of the gear transmission system at home and abroad mostly focused on the spur gear and spur planetary gear system. In the reports on the nonlinear dynamics of double-helical gear, there are few reports on the nonlinear frequency response characteristics of double-helical gear. Based on the above research, a torsional dynamic model of the doublehelical gear pair system is established in this study. e timevarying meshing stiffness, static backlash, dynamic backlash, static transmission error, external dynamic excitation, and other factors are considered in the model. e influence of parameter changes on the nonlinear dynamic characteristics of the double-helical gear pair system is studied by using the incremental harmonic balance method. e research results can provide reference for vibration reduction analysis and structure optimization of double-helical gear.

Torsional Nonlinear Dynamic Model of Double-Helical Gear Transmission
As shown in Figure 1, the torsional nonlinear dynamic model of double-helical gear pair is shown. e time-varying meshing stiffness, viscoelastic damping, backlash, static transmission error, and external dynamic excitation are considered in the model. According to Newton's second law, the torsional nonlinear dynamic equation of double-helical gear pair can be written as In the above formula, I ij (i � p, g.j � L, R) represents the polar moment of inertia on the equivalent nodes O pL , O pR , O gL , and O gR at the center of each helical gear, θ ij (i � p,g; j � L, R) represents the angular rotation displacement of the equivalent node at the center of each helical gear, r ij (i � p, g; j � L, R) represents the radius of the base circle of each helical gear, β j (j � L, R) represents the helix angle of the left and right helical gear pairs, k mj , c mj (j � L, R) represents the meshing stiffness and damping of the left and right helical gear pairs, δ j (j � L, R) represents the relative meshing displacement of the left and right helical gear pairs, e j (j � L, R) represents the static transmission error of the left and right helical gear pairs, k it ,c it (i � p, g) represents the torsional stiffness and torsional damping of the intermediate shaft of double-helical gear, and T pL (t), T pR (t), T gL (t), T gR (t) represent the input torque and output torque acting on the driving helical gear and the driven helical gear, respectively.
Each moment can be expressed as follows: 2 Shock and Vibration where T pL , T pR , T gL , and T gR are the external nominal moments, and T pLa (t), T pRa (t), T gLa (t), and T gRa (t) are the external dynamic moments. e four degrees of freedom θ pL , θ pR , θ gL , and θ gR in the nonlinear torsional dynamic equation of double-helical gear transmission are not independent of each other. e relative displacement method is used to eliminate the rigid body displacement of the system to solve the differential equation. e absolute torsional angular displacement θ pL , θ pR , θ gL , and θ gR on the original mass nodes is replaced by the relative meshing displacement δ L , δ R along the meshing line direction of the left and right helical gears and the relative torsional line displacement c P , c G between the coaxial adjacent helical gears. e relative meshing displacement of the left and right helical gear pairs along the direction of the meshing line is expressed as follows: e expression of relative torsional line displacement of adjacent coaxial helical gears is as follows: e actual meshing displacement of double-helical gear along the meshing line is expressed as a piecewise function as follows: where 2b is the length of backlash. For the actual doublehelical gear pair transmission, with the difference between the relative meshing displacement along the meshing line direction and the backlash, that is, the actual meshing displacement along the meshing line direction, the doublehelical gear pair transmission may be in three different working states, namely, no impact, single-side impact, and double-side impact [22,23,27], as shown in Figure 2.
In the existing research, the backlash function is expressed as piecewise linear. Generally, the backlash of gear teeth is obtained by static measurement, and the change of backlash in the process of gear impact is not considered, so the simulation results cannot truly reflect the transmission characteristics of the gear system. In order to reveal the influence of the actual tooth surface characteristics on the vibration and noise of the gear transmission system, the constant backlash model and dynamic backlash model will be considered simultaneously in this study. e constant backlash model takes the backlash as a constant according to the traditional modeling method. e expression is as follows: e dynamic backlash model assumes that the backlash is in the form of harmonics. e expression is as follows: where b 0 is the equivalent backlash, W is the backlash variation describing the wear of gear tooth, and b l is the amplitude of dynamic backlash representing the deformation of the supporting bearing and axle. By substituting  (9) into (2), the dynamic equation (8) is obtained: where m c is the equivalent mass of the left and right helical gear pairs, and f a is the fluctuation of external load excitation, which are, respectively, expressed as follows: where f l is the fluctuation coefficient. e time-varying meshing stiffness function and static transmission error function are expressed in the form of Fourier series [2-6, 8, 28, 29]: where k mean is the average stiffness, ε l is the fluctuation coefficient reflecting the change of stiffness, and e l is the fluctuation coefficient reflecting the change of static transmission error. When the excitation of static transmission error is maximum, the time-varying meshing stiffness should be the minimum, that is, the two have inverse phase relationship. e time-varying meshing stiffness curve is shown in Figure 3. e meshing damping c mj and torsional damping c it are calculated by the following formula: where ζ is the damping ratio of the system. According to the knowledge of material mechanics, the calculation formula of torsional stiffness of intermediate shaft section is as follows: where G is the shear elastic modulus of the gear material, J is the polar moment of inertia of the cross-section, l is the length of the intermediate shaft segment, and D and d are the outer diameter and inner diameter of the shaft, respectively. Equation , ω n is the natural frequency, t is the normalized time, and Ω is the normalized meshing frequency. In this study, the derivative of time t in equations (8) is transformed into the derivative of dimensionless time t, and finally, the dimensional unified dynamic equation is obtained: Shock and Vibration 5 e stiffness term in the above formula is as follows: e static transmission error term in the excitation vector is as follows: e dynamic excitation of external torque is as follows: For the piecewise backlash function, if b c � b 0 is taken in the dimensionless process, the two kinds of backlash models are transformed into Equation (20) is written in the matrix form as follows: where M is the dimensionless mass matrix, C is the dimensionless damping matrix, K is the dimensionless stiffness matrix, F(t) is the dimensionless excitation vector, and q(t) is the dimensionless displacement vector. Limited to space, only the stiffness matrix is listed as follows:

Approximate Solution Based on the Incremental Harmonic Balance Method
By letting a new time scale τ � Ωt, equation (22) can be represented as e Newton-Raphson algorithm is used to solve the piecewise linear difference system represented by equation (24), and the solution of the differential equation is obtained as follows: where q 0 (τ) is the approximate solution of equation (24), and Δq(τ) is the incremental equation. Let any term q j0 in the n-order approximate solution of formula (25) and any term Δq j in the incremental equation are, respectively, Δa jn cos(nθ) + Δb jn sin(nθ), where a jn , b jn , Δa jn , Δb jn are the Fourier coefficients, and N is the number of terms of the Fourier series. e nonlinear piecewise linear function of backlash may be expressed by a first order Taylor expansion as By substituting equations (26) and equation (29) into equation (24), Taylor series expansion is carried out, and the high-order term is omitted, and equation (24) becomes linearized as (lΩ) 2 e l sin(lτ) Order j0 , a j1 , . . . , a jN , b j1 , . . . , b jN T , By using the Galerkin process [5,8,[30][31][32], the left and right ends of equation (29) are multiplied by cosiτ, sin(iτ) (i � 0, 1, 2, 3, . . ., N). e linear equations are obtained by integrating between 0 and 2π.
Equation (36) is the iterative calculation formula of the steady-state periodic response of the nonlinear system with ΔA as unknown quantity, which is derived by using the incremental harmonic balance method.
According to equation (36), the approximate solution with arbitrary precision can be obtained, that is, given the initial value of iteration, ΔA is obtained; thus, the initial value A of the next iteration can be obtained. When the accuracy requirement is met, the iteration will be terminated.

Comparison between Analytical Solution and Numerical
Solution. e design parameters of double-helical gear pair are given in Table 1.
In order to fully study the nonlinear characteristics of double-helical gear under the coupling effect of time-varying meshing stiffness and backlash, the amplitude frequency characteristic curves of the system under four different models are obtained, and the results are shown in Figure 4.
In Figure 4, the ordinate represents the maximum Δ Lmax and minimum Δ Lmin nonlinear vibration displacement of the double-helical gear system along the meshing line, and the abscissa represents the dimensionalized frequency Ω. Four different models are, respectively, the constant backlash model under combined internal and external dynamic excitations as shown in Figure 4(a), the constant backlash model under external dynamic excitation as shown in Figure 4(b), the constant backlash model under internal dynamic excitation as shown in Figure 4(c), and the dynamic backlash model under combined internal and external dynamic excitations as shown in Figure 4(d). e solid line represents the result obtained by incremental harmonic balance method, and the circle represents the result obtained by the Runge-Kutta method with variable step size of fourth and fifth order. Four different models are as follows: Figure 4 shows that there are multivalue solutions and jumping phenomena in the double-helical gear system under four different models. e multivalue solution interval of the constant backlash model is 0.78-0.86 and that of the dynamic backlash model is 0.83-0.84; the system includes no impact state and single-side impact state, no double-side impact state is found, and tooth back separation does not occur. In the high frequency range of dimensionless frequency from 1.794 to 1.803, transient instability occurs in the numerical method but not in the analytical method. At the same time, it can be found that the nonlinear dynamic characteristics under combined internal and external excitations are not simply obtained by linear superposition of internal and external excitations, respectively. e multivalued properties and jump phenomena are corresponding with no impact motion and single-side impact motion. Figure 5 shows the steady-state response of the gear pair system when Ω � 0.83. e basic system parameters are L � 3, N � 11, F 1 � 1.3534, ζ � 0.045, ε 1 � 0.16, ε 2 � 0.07, ε 3 � 0.03, e 1 � 0.03, e 2 � 0.01, e 3 � 0.01, f 1 � 0.06, f 2 � 0.03, f 3 � 0.01, and W � 0. e phase diagram and time history diagram of main resonance response and superharmonic response of the doublehelical gear system are obtained by using the incremental harmonic balance method. Meanwhile, the Runge-Kutta numerical method is used for comparison and verification. e results are shown in Figure 6. Figure 6 shows that when Ω � 1, the phase diagram is approximate ellipse, the system is a single period harmonic response, and the gear is in the state of no impact; when Ω � 1/2 and Ω � 1/3, the phase diagram is a noncircular closed curve, the system response is superharmonic response, and the meshing state of gear teeth is not changed. e generation of the superharmonic response solution of the system is caused by the high-order harmonic component   Figure 7.
From the comparison between the linear time-varying model (zero backlash) and the nonlinear time-varying model (with backlash) in Figure 7, it can be seen that the existence of backlash nonlinearity leads to the nonlinear dynamic characteristics such as multivalued solution and amplitude jump of the double-helical gear pair system. e basic system parameters are L � 3, N � 11, F 1 � 1.3534, ζ � 0.045, ε 1 � 0.16, ε 2 � 0.07, ε 3 � 0.03, e 1 � 0.03, e 2 � 0.01, e 3 � 0.01, f 1 � 0.06, f 2 � 0.03, and f 3 � 0.01. e influence of different constant backlash and dynamic backlash amplitude on the amplitude frequency characteristics of the system is shown in Figure 8.
For constant backlash, five models of constant backlash are considered: (1) W � − 0.6; (2) W � − 0.3; (3) W � 0; (4) W � 0.3; and (5) W � 0.6. It can be seen from Figure 8(a), with the increase of the constant backlash, the vibration displacement amplitude of the gear system increases, the multivalue solution region is wider, and the jumping phenomenon still exists. For dynamic backlash, four kinds of dynamic backlash models are considered: It can be seen from Figure 8(b) that with the increase of the dynamic backlash amplitude, the vibration displacement amplitude of the gear system slightly increases, the jumping phenomenon still exists, and the multivalue solution interval gradually decreases until it disappears. It can be seen that increasing the dynamic backlash can control the collision behavior of the gear system.

Influence of Time-Varying Meshing Stiffness on Amplitude
Frequency Curve of the System. e basic system parameters are L � 3, N � 11, F 1 � 1.3534, ζ � 0.045, e 1 � 0.03, e 2 � 0.01, e influence of time-varying meshing stiffness on the amplitude frequency curve of the system under the constant backlash model and dynamic backlash model is shown in Figure 9.
In addition, compared with the constant backlash model, the dynamic backlash model can control the nonlinear vibration of the system more effectively by reducing the fluctuation coefficient of the meshing stiffness.  Figure 10.

Influence of Damping on
It can be seen from Figure 10 that damping can restrain the step of amplitude frequency characteristic curve, the emergence of multivalue solution interval, and the amplitude of nonlinear vibration displacement. Increasing the system damping can effectively control the nonlinear dynamic response of the gear system; compared with the constant backlash model, increasing the damping ratio under the dynamic backlash model can more effectively control the nonlinear vibration of the system.
is study considers three kinds of static transmission error amplitudes: (1) e 1 � 0.1, e 2 � 0.05, e 3 � 0.03; (2) e 1 � 0.03, e 2 � 0.01, e 3 � 0.01; and (3) e 1 � 0, e 2 � 0, e 3 � 0. It can be seen from Figure 11 that the amplitude frequency characteristic curve of the system changes greatly with the change of the amplitude of the static transmission error, when e 1 � 0, e 2 � 0, and e 3 � 0, and the amplitude frequency curve still has a jump phenomenon and multivalue solution phenomenon, which indicates that the static transmission error in doublehelical gear is not the main reason for the nonlinearity of the system. When e 1 � 0.1, e 2 � 0.05, e 3 � 0.03 and e 1 � 0.03, e 2 � 0.01, e 3 � 0.01, the larger the amplitude frequency curve of the system appears, and the larger the value of static transmission error amplitude is, the greater the amplitude frequency characteristic curve amplitude is.

Influence of Excitation Force Amplitude on Amplitude
Frequency Curve of the System. e basic system parameters are L � 3, N � 11, ζ � 0.045, ε 1 � 0.16, ε 2 � 0.07, ε 3 � 0.03, e 1 � 0.03, e 2 � 0.01, e 3 � 0.01, f 1 � 0.06, f 2 � 0.03, f 3 � 0.01, W � 0, b 1 � 0.05, b 2 � 0.02, and b 3 � 0.01. In the constant backlash model and dynamic backlash model, the influence of excitation amplitude force on the amplitude frequency curve of the system is shown in Figure 12. It can be seen from Figure 12 that with the increase of the excitation force amplitude, the amplitude frequency characteristic curve amplitude of the system increases, the jump phenomenon still exists, and the multivalue solution phenomenon disappears slowly. When F 1 � 3.3534, the amplitude frequency characteristic curve amplitude is the largest, while there is no multivalue solution interval. When F 1 � 1.3534, the amplitude of amplitude frequency characteristic curve is the smallest, and the range of unstable solution interval is the largest. e transition frequency of curve step basically does not change, and the results of the two backlash models have little difference. With the increase of the excitation amplitude, the amplitude frequency characteristic curve of the system still contains the no impact state and the single-side impact state, and the increase of the excitation amplitude has little effect on the change of the meshing impact state of the double-helical gear.

Conclusions
In this study, the torsional dynamic model of double-helical gear transmission is established. e factors such as timevarying meshing stiffness, static transmission error, constant backlash, dynamic backlash, and external dynamic excitation are considered in the model. e frequency response characteristic curve of the system is better solved by using the incremental harmonic balance method. It is found that the system has multivalue solution and jumping nonlinear characteristics, and the results are in good agreement with the numerical integration results. e conclusions are as follows: (1) ere are main resonance response and superharmonic response in the double-helical gear system. e generation of the superharmonic response solution is caused by the high-order harmonic component in time-varying meshing stiffness and dynamic excitation. backlash amplitude, the vibration displacement amplitude of the gear system increases slightly, and the multivalue solution range decreases or even disappears.
(3) e results show that the time-varying meshing stiffness and static transmission error can stimulate the amplitude frequency curve of the system, and the damping can restrain the amplitude frequency curve of the system. Compared with the constant backlash model, the vibration of the double-helical gear system can be further controlled by changing the system parameters in the dynamic backlash model.

Data Availability
e data used to support the findings of this study are included within the article.

Conflicts of Interest
e authors declare that they have no conflicts of interest.