Abstract
Bevel gears are widely used in aerospace transmission systems as well as modern mechanical equipment. In order to meet the needs and development of aerospace, high-speed dynamic vehicles, and various defense special equipment, higher and higher requirements are made for the high precision and stability of gear transmission systems, as well as the prediction and control of noise and vibration. Considering the nonlinear factors such as comprehensive gear error and tooth side clearance, a dynamic model of the three-stage gear transmission system is established. The relevant physical parameters, geometric parameters, and load parameters in the gear system are considered random variables to obtain the stochastic vibration model. When the random part of the random parameters is much smaller than the deterministic part, the vibration differential equation is expanded into a first-order term at the mean of the random parameter vector according to the Taylor series expansion theorem, and the ordering equation is solved numerically. Based on the improved stochastic regression method, the nonlinear dynamic response analysis of the three-stage gear train is carried out. This results in a relatively stable system when the dimensionless excitation frequency is in the range of 0.716 to 0.86 and the magnitude of the dimensionless integral meshing error is < 1.089.
1. Introduction
Gearing systems are widely used in various types of machinery and equipment, and ensuring the smoothness and accuracy of the transmission process has always been the goal of industrial production [1]. With the development of wind turbine drive train power towards the megawatt level, the accuracy of the dynamic response analysis of the three-stage gearbox drive train model has been increasing. Therefore, the study of the coupled nonlinear dynamic characteristics of large generator gearboxes under variable loads, and then provide important theoretical value and practical engineering significance for the design of three-stage gearboxes [2]. In recent years, many domestic and foreign scholars have conducted in-depth studies on the performance research methods and dynamics models of three-stage gearbox drive systems, and have made certain achievements [3, 4, 5]. The coupled nonlinear dynamics model of the gear-rotor-bearing-box of a large three-stage gearbox with multi-stage planetary gearing was established under the combined effect of internal and external load excitation, and the vibration response analysis of the three-stage gearbox was carried out by using the modal superposition method. Reference [6] accurately modeled the gearbox drive train of a large wind turbine and studied its drive train under variable wind loads. Reference [7] established a bending-torsion coupled dynamics model of a megawatt-class wind turbine gearbox drive train and calculated the structural sketch of the three-stage gearbox drive train using the fourth-order Runge–Kutta method.
Wind turbine gearboxes can be configured in different types, such as three-stage fixed shaft gear systems, planetary gear systems, and planetary and parallel shaft gear systems. The form of transmission mechanism used in the three-stage gearbox depends on the capacity requirements of the three-stage unit. The development of offshore three-stage technology and the need to reduce the cost of three-stage have strongly promoted the development of large single-capacity three-stage units. Three-stage planetary gearboxes are widely used in high-power double-fed three-stage gearboxes because of their strong load-bearing capacity, large speed increase ratio, and compact structure. Compared with single-stage planetary gearboxes, three-stage planetary gearboxes are more susceptible to vibration and noise and failure due to interstage coupling and three-stage transmission of torque [8]. Therefore, it is important to analyze the dynamic characteristics of three-stage planetary gear trains. A series of research results have been obtained for the analysis of the dynamics of a three-stage planetary gear train in three-stage gearboxes.
The dynamics of a single pair of a gear pair system has been studied more systematically. In [8], a time-varying nonlinear dynamics model was established, and the effect of side clearance on the meshing state and its frequency response characteristics was considered [9]. The parameters of the single-stage spur gear system are studied in detail, and the effects of parameter variations on the dynamic characteristics are presented. In [10], an approximate analytical solution of the three-stage nonlinear gear train was carried out using the harmonic balance method, and the effect of the variation of key parameters on the system response was analyzed [11]. The dynamic simulation technique of a three-stage gear train considering only the time-varying meshing stiffness was studied [12]. The nonlinearities of the intermediate shaft system and the idler system with time-varying meshing stiffness are investigated, and the analytical and numerical solutions are compared. In [13], the dynamic characteristics of the three-stage idler system with nonlinear parameter excitation are presented, and the periodic steady-state solutions are obtained by both analytical and numerical methods.
Tooth surface errors are unavoidable in the gear manufacturing process. In this paper, a modified stochastic regression method is used to solve the system of vibration differential equations for a three-stage gear train, and the effects of the changes of the integrated dimensionless excitation frequency and the integrated meshing error on the dynamic response of the system are analyzed. It is found that the system is in a relatively stable state when the dimensionless excitation frequency is in the range of 0.716–0.86 and the magnitude of the integrated dimensionless meshing error is <1.089.
2. Related Work
In order to meet the practical needs of low speed and a heavy load and a large transmission ratio, the study of the dynamics of the gear transmission system is becoming a hot spot for engineering applications and academic research. According to the research focus, [14] used the finite element method, the centralized parameter method, and the experimental modal analysis method, respectively. In [15], the finite element model and the centralized parametric model of the gear set were established, and the results of the two models were compared and analyzed, and the results of the two models were found to be consistent. In [16], the dynamic characteristics of the system under different fixation methods of the solar wheel were analyzed. In [17], the vibration model of the planetary gear was developed by the centralized parameter method considering the meshing phase relationship, and the modalities and vibration modes of the system were calculated for different modalities. In terms of nonlinear dynamics of planetary gear sets, [18] analyzed the dynamics of planetary gears under changing meshing stiffness by using rectangular waves to represent the time-varying meshing stiffness of the gear teeth. In [19], the effect of manufacturing errors on the nonlinear dynamics of planetary gears were studied, and the differences in the dynamics of the system under normal conditions, in the presence of tooth profile errors, and when centrifugal forces were considered and compared. In [20], the dynamics of planetary gears under strong nonlinear factors with multiple clearances is considered, and it is important to note that the analytical solution of the system vibration differential equations is obtained by the harmonic balance method, which is a typical application of the approximate analytical solution method under strong nonlinear factors. The study of the concentrated mass model of planetary gear sets focuses on the pure torsional vibration model and the meshing coupling model. Theoretical tooth surface and real tooth surface STE are introduced into the dynamics model in the form of discrete points, and the effects of real tooth surface STE and theoretical tooth surface STE on the dynamics performance of the curved bevel gear pair are compared and analyzed.
3. Complex Modal Structural Systems
The basic equation for the vibration of a gear system and the complex modal forced vibration equation for an N-dimensional linear system can be expressed aswhere M, C, and K denote the N ×N order mass matrix, damping matrix, and stiffness matrix, respectively, X and Q(t) denote the generalized displacement vector and external excitation, respectively, denotes the generalized velocity vector and denotes the generalized acceleration vector. The free vibration equation of the complex modal structural system can be expressed as
It may be assumed that, by substituting into (2), the corresponding right eigenvalue problem is obtained as follows:
The corresponding concomitant eigenvalue problem is , and the left eigenvalue equation is obtained by transposing the concomitant eigenvalue equation as follows:where the vectors x and y denote the right and left eigenvalues, respectively, and λ denotes the eigenvalue of the vibration equation.
To simplify the above representation, a right state vector u is introduced to replace , satisfying , where T is the state transition matrix and can be expressed as . Similarly, another left state vector is used instead of , satisfying . Thus, (3) and (4) can be rewritten as
Among them
After the state transformation, (6) can be rewritten as
(7) and (8) have the same eigenvalues and can be derived from the following equation:
(8) is an algebraic equation with 2N eigenvalues in the complex domain, which can be expressed as .
A and B can be expressed as
4. Matrix Ingestion Method for Complex Modal Structural Vibrations
Since the non-diagonal of the consistent mass matrix contains non-zero elements, when vectorial finite elements are used for dynamic analysis, the consistent mass matrix is often replaced by the concentrated mass matrix in order to reduce computational effort and storage space. The most common method of constructing a centralized mass matrix is to distribute the mass of a cell equally among the connected nodes. When the mass of the unit or the internal nodes are not uniformly distributed, a higher-order difference function is used for mass distribution.
For the rod system structure, the elements on the main diagonal of each row of the centralized mass matrix element on each row of the consistent mass matrix is the sum of all elements in the corresponding row of the consistent mass matrix, i.e.,where is the concentrated mass matrix. is the consistent mass matrix. is the th-order square matrix composed of 1. Substituting equation (5) into (6) further simplification yields
In practical engineering, the variation of structural parameters is reflected in the variation of the mass matrix. Considering the small perturbations of the parameters around their mean values, the mass matrix can be expressed as
According to the matrix regression method, the eigenvalues and eigenvectors of a complex modal structural vibration system can be expanded into the form of a regression series:
Substituting (11) and (13) into (12), we get
Similarly, substituting equations (12) and (13) into (14) gives
Expanding (15) and ignoring the terms of O(), comparing the coefficients of the same power of ε on the left and right sides of the equation yieldswhere ε is a small parameter and has ε = 0 for the original system; denotes the mean values; and denotes the first-order regression term of , respectively. Further, the same operation of (13) and (14) for (17) giveswhere and represent the first-order and second-order registers of the th eigenvalue, respectively. In addition, the first-order and second-order registers of the left eigenvector and the first-order and second-order registers of the right eigenvector can also be found by (15).
Expanding on the original right state eigenvector yields
Substituting (17) into (18) yields
The left multiplication of (19) by has
Considering the modal orthogonal normalization conditions expressed in (14) and (15), (20) can be reduced to
Only if , then can be obtained from (22)
(20) can be further rewritten as
5. Uncertainty Propagation Analysis of Random Complex Eigenvalues
The mass matrix, damping matrix, and stiffness matrix of real engineering structures are generally considered as stochastic structural parameters due to the variation of raw material properties of gears, manufacturing errors, and loading environment. In this paper, the matrix regression method is introduced to the eigenvalue analysis of complex modal structures with stochastic parameters, and also provides a basis for uncertainty propagation analysis of asymmetric systems. Considering the stochastic nature of the structural parameters, the stiffness matrix K, the mass matrix M, the damping matrix C, the complex eigenvalues , the left complex eigenvector y, and the right complex eigenvector x can be expressed in the following form:where ε is a small parameter, , , and are the deterministic and random components of K, M, C, and , respectively. Also, it is assumed that the mean value of is zero.
6. Numerical Solution of the Model
In this paper, we solve the differential equations of the system in terms of the magnitudes and vibrations, and briefly compare the results of the two numerical integration methods. Given the initial values of the system vibration, the first 800 transient response periods are omitted, and the dynamic response of the system is programmed in MATLAB to analyze the effect of the variation of excitation frequency, mesh damping ratio, and tooth side clearance on the bifurcation characteristics of the system. The basic parameters of the two-stage planetary gear train of the three-stage gearbox are as follows: rated power P = 2500 kW, normal speed range 16r/min∼21r/min, transmission ratio = 5.25, = 5.28, number of teeth of the first stage Zs = 31, Zpi = 47, Zr = 125, number of planetary wheels = = 3, and the basic parameters of each gear are shown in Table 1.
The presence of static transmission error can aggravate the meshing shock of the gear teeth and cause abnormal vibration of the system, therefore, this paper discusses the effect of the integrated meshing error on the bifurcation characteristics of the system.
The kinetic analysis of the solution results by the global bifurcation diagram and maximum Lyapunov exponent (LLE) is shown in Figures 1 and 2. Given the dimensionless excitation frequency = 0.86, the mesh damping ratio ξ = 0.07 and the tooth side clearance b = 1, the effect of the integrated meshing error on the dynamic characteristics of the system is discussed with the dimensionless integrated meshing error magnitude E as the bifurcation parameter (for simplicity, the image coordinates are denoted as and as E). According to Figure 1, the three-stage planetary gear system has a rich nonlinear dynamic behavior when the magnitude of the integrated dimensionless meshing error is used as the bifurcation parameter. Combined with Figure 2, it is found that the system is relatively stable at low values of the integrated meshing error, when LLE is negative, and it is clear that the system enters the chaotic response as a multiperiodic bifurcation.


The dynamics of the two-stage planetary gear and the one-stage parallel shaft gear system are further verified and analyzed by combining the phase trajectory diagram and the Poncelet cross section diagram with the aid of the time history diagram and the spectrum diagram, as shown in Figure 3. When the dimensionless integrated meshing error value = 0.05, the system is in a single-cycle motion state, the time history graph is close to a sinusoidal curve, the vibration amplitude between cycles is basically the same, the phase trajectory line is an elliptical closed curve, the Poincaré diagram is a single point and the fundamental frequency signal in the spectrum is very prominent at this time, indicating that the system is in a typical single-cycle state.

With the increase of the integrated meshing error, when 0.779, as shown in Figure 4, the time course diagram changes periodically with two amplitudes as one cycle, the phase trajectory line is concave from one elliptic curve to two elliptic curves, the Poincaré cross-sectional diagram splits from one single point to two single points, and the FFT spectrum is shown as octave components, indicating that the system is in a two-cycle motion state. Combining Figures 3 and4, it can be found that the motion state of the system splits and changes from single-cycle to two-cycle and four-cycle with the increase of the magnitude of the integrated dimensionless mesh error. At the integrated dimensionless meshing error (1) 1.259, strange attractors appear in the cross-sectional plot of Poincaré, indicating that the system enters the chaotic state through multi-cycle bifurcation at this time. When (1) reaches 0.939, the system bifurcates from the two-cycle state to the two-cycle state with an LLE value of -0.001713. Therefore, the “negative pole” that appears in the plot of the maximum Lyapunov exponent of the system under the change of the integrated mesh error can characterize the dynamics of the two-cycle bifurcation, as shown in Figure 5.


7. Conclusion
In this paper, a modified stochastic regression method is used to solve a system of vibration differential equations for a three-stage gear train, and the effects of changes in the integrated dimensionless excitation frequency and the integrated meshing error on the dynamic response of the system are analyzed. A new static no-load transmission error model is established, which has its unique advantages in the fault diagnosis of gears. The model no longer requires tooth surface equations, avoiding the difficulty of merging tooth surface errors into tooth surface equations. The effect of the static unloaded transmission error of the real tooth surface on the dynamic characteristics of the bevel gear is analyzed. The model is also applicable to other types of gears. The global bifurcation diagram, Poncelet cross section diagram, and phase trajectory diagram are used as the main basis, and the time history diagram and the spectrum diagram are used as auxiliary tools to determine the range of bifurcation parameters when the three-stage gearing system is in a relatively stable state. It is found that the system is in a relatively stable state when the dimensionless excitation frequency is in the range of 0.716–0.86 and the magnitude of the dimensionless integrated meshing error is < 1.089.
Data Availability
The data used to support the findings of this study are available from the corresponding author upon request.
Conflicts of Interest
The authors declare that they have no conflicts of interest regarding this work.