# Research and Application of Data-Driven Modeling Method for Gear System Based on Time-Varying Characteristics

Authors: Zi-wen Zhu; Xiaojun Zhu; Yong-bo Pan
Journal: Mathematical Problems in Engineering
 (2022)

---

## Abstract

The time-varying characteristics of the gear system have an essential influence on its vibration and stability characteristics. Aiming at this characteristic and taking the generalized force caused by load torque as input and dynamic transmission error as output, a data-driven modeling method for gear time-varying system is studied using the periodic time-varying system identification theory. The effectiveness of the proposed method is verified by the lumped mass model of the gear system. This method has high modeling accuracy and can accurately characterize the time-varying characteristics of the gear system. Then, the virtual experimental platform of gear system dynamics is established based on the finite element theory. The method successfully extends to the virtual experimental platform, laying a foundation for analyzing the gear system’s dynamic characteristics.


---

## Body

## 1. Introduction

Gears are widely used in manufacturing. Gear systems are complex, elastic structural systems, and the study of their dynamic characteristics has been of great interest. However, the dynamic modeling of the gearbox is recognized as one of the challenging topics. From the perspective of system identification theory, it is essential to research the modeling of the dynamic characteristics of gear transmission systems.The dynamic characteristics or parameter identification of gearboxes has greatly interested scholars. Tong et al. [1]proposed a multi-input single-output (MISO) model to identify the local resonance excited by gear meshing impact and verified the method’s performance through experiments. Dai et al. [2] evaluated the meshing force. Based on the extended frequency response function (FRF) concept, the planetary meshing force identification model was established. Amthor et al. [3] established the nonlinear analytical model of the gear system, deduced the parameter identification algorithm of the model based on theoretical modeling, and verified the method through experiments. Sawalhi and Randall [4] identified the number of teeth and gear speed by measuring the gear vibration acceleration signals. Pedersen et al. [5] applied the modal analysis method to analyze the gear system with periodic time-varying (PTV) meshing stiffness and obtained the fundamental frequency and resonant frequency variation laws of the gear with the rotational speed. Based on the modal analysis theory, Saxena et al. [6] studied the influence of meshing stiffness of healthy gear and cracked gear on the gear system’s natural frequency, modal shape, and FRF. Ericson and Parker [7] used the classical experimental modal analysis technique to characterize the dynamic characteristics of planetary gear, and the natural frequency, vibration mode, and FRF of planetary gear are obtained. From the measured FRF, Aykan et al. [8] identified the gear system nonlinearity and obtained the corresponding mathematical model. Mbarek et al. [9] studied the influence of different load conditions and meshing stiffness fluctuations on planetary gears’ modal and natural frequency. Dong et al. [10] analyzed the frequency response characteristics of the gear system under constant backlash and dynamic backlash by incremental harmonic balance method and further verified the results by a numerical method. Yang et al. [11] studied the FRF of a nonlinear time-varying cylindrical spur gear system under multiple excitations, which provides a valuable reference for reducing the vibration of the gear system. Wang et al. [12] proposed a parameter identification method based on wavelet transform and correlation filtering and applied this method to the gearbox. The time-varying characteristics of gear systems have not been systematically studied.Gear (dynamic) transmission error is one of the leading indicators to characterize the dynamic properties of a gear system. Numerous scholars have done much work on the simulation calculation and experimental measurement of transmission errors. Wang et al. [13] established a finite element model to predict the dynamic transmission error (DTE) by considering gear eccentricity and load change. Xiang and Gao [14] analyzed the mechanism of manufacturing and assembly errors, proposed a transmission error analysis method for helical gears with errors, and verified the method by experiments. Kong et al. [15] analyzed the influence of gearbox shell flexibility on the DTE of the gear transmission system. Duan et al. [16] systematically studied the transmission error of the gearbox by establishing the rigid-flexible coupling dynamic model of the gearbox and establishing a transmission error test bench. Wei et al. [17] analyzed the influence of contact ratio, support stiffness, meshing damping, the backlash on DTE, and the internal relationship between DTE and the above parameters. Anichowski et al. [18, 19] integrated the DTE measurement system with an accelerometer into the high-speed gear tester. They measured the DTE of gear pairs with different indexing error forms. Benatar et al. [20] measured the transmission errors of helical gears with different profiles at low speeds (quasi-static) and dynamic conditions based on the encoder and acceleration sensor. The experimental results showed that the tooth modifications significantly affected the amplitude of the transmission system’s static and dynamic transmission errors. Xiong and Chen [21] investigated the effect of modification and misalignment gears on DTE under different load and drive speed conditions using high-precision optical encoders. In addition, the relationship between DTE and manufacturing errors had also been studied in depth [22–25].The paper is structured as shown in Figure1. First, taking the torque generalized force transmitted by the gearbox as the input and the DTE as the response, the gear system’s time-varying response function identification method is established by using the PTV identification theory. Then, the effectiveness of the proposed method is verified by the time-varying dynamic numerical model of the gear system [26]. Finally, the method is extended to the virtual experimental platform of the gear system to identify its frequency response function. This paper presents a new data-driven modeling method for dynamic research of gear systems, which can provide a basis for vibration and noise reduction, fault diagnosis, and intelligent control of gear transmission systems.Figure 1
Main structure of this paper.


## 2. Identification Method of a Periodic Time-Varying System

### 2.1. Periodic Time-Varying System

PTV systems are a class of common time-varying systems in practical engineering, such as electric motors, fans, shafts, and helicopter blades. The dynamic properties of these systems often contain a PTV coefficient.Since the number of teeth involved in meshing simultaneously in the gear transmission system changes periodically with time, the gear’s meshing stiffness, meshing damping, and other parameters will change periodically with time [27]. Therefore, a gear system can be considered a PTV system [28].Under the condition of uniform speed, the modeling and identification of the time-varying FRF of the gear system are carried out using the PTV system identification theory [29]with the torque generalized force and DTE as the input and output data of the system, respectively.


### 2.2. Input-Output Model of a Periodic Time-Varying System

The linear PTV systemG can be well described by a parallel structure consisting of infinitely many weighted linear time-invariant systems, as shown in Figure 2(a). In the figure, the triangle blocks denote a time-domain multiplication. Moreover, the (triangles) gains are written in a shorthand notation as lωsys instead of ejlωsyst. In practice, we use the branch structure composed of finite (such as Nb)weighted linear time-invariant systems to better approximate system G(1)Gω,t=∑l=−Nb+NbGlωejlωsyst,where Gω,t is called the instantaneous transfer function (ITF), Glω is called harmonic transfer function (HTF), and ωsys is the time-varying frequency of the system, namely, the meshing frequency of the gear transmission system.Figure 2
Principle diagram of a linear periodic time-varying system. (a) System direct equivalent model. (b) System indirect equivalent model.
(a)(b)In this paper, the ITF is identified in the frequency domain. The following are the basic assumptions needed in the identification.Without changing the input-output relationship of the system, we can exchange the time-varying gain module and the linear time-invariant module in Figure2(a) to obtain the model as shown in Figure 2(b). This model can be regarded as a linear time-invariant system with multiple inputs and a single output, and the input-output relationship in the frequency domain of the model can be expressed as follows.(2)Y0k=∑l=−Nb+NbHlωkU0k−pl,where U0k, Y0k are the spectra of discrete Fourier transform of input and output vectors. Hlωk is the HTF corresponding to the model in Figure 2(b), p=ωsys/ω0, and ωk=kω0.


### 2.3. Local Polynomial Estimation of the ITF

For the abovementioned model, this paper uses the local polynomial method (LPM). The LPM is a nonparametric identification method commonly used for linear time-invariant systems with multiple inputs and multiple outputs. This method uses polynomials to model FRF and transient terms in the local interval and uses the least square method to estimate model parameters based on smooth frequency characteristics [30].Assuming that the noise model is additive, according to (2), we can express the input-output relationship of a linear PTV system with output noise as follows.(3)Yk=HωkUk+WkUk=U0k.Uk is the input vector,(4)Uk=Uk−pNb…Uk…Uk+pNbT.The FRF between inputU0k and output Y0k is as follows.(5)Hωk=HNbωk…H0ωk…H−Nbωk.In addition,Wk is an error term consisting of the output noise (NY) and the output stochastic nonlinear disturbance (YS).In a local frequency band, that is, at the frequency linesk+l for (l=−n,−n+1,…,0,…,n), the FRF Hωk can be locally approximated by a low-order polynomial.(6)Hωk+l=Hωk+∑r=1Rhrklr+Ol/NR+1,where R is the order of the polynomial.Neglecting the bias term in (6), at the frequency line k+l, (3) can be expressed as follows.(7)Yk+l=ΘkKk+l+Vk+l,where(8)Θk=Hωk,h1k,h2k,…,hRk,is a 1×R+12Nb+1 row vector of the unknown parameters, and the vector Kk+l is defined as follows:(9)Kk+l=1l...lR⊗UTkT,where ⊗ is the Kronecker product operator. Evaluating (7) at l=−n,−n+1,…,0,…,n and grouping the equations in a matrix equation finally yield as follows:(10)Yk=ΘkKk+Vk,where Yk, Kk, and Vk are, respectively, 1×2n+1, R+12Nb+1×2n+1 and 1×2n+1 matrices.Through the least square method, we yield as follows.(11)H^ωk=YkSk,(12)S=KHKKH−1I2Nb+1,0H,where frequency argument k in (12) is removed for notational simplicity and xH is the complex conjugate transpose of x.The estimation of the HTFG^lω can be obtained from equations (5), (11), (13)(13)G^lωk=H^lωk+pl,where p=ωsys/ω0.Finally, the nonparametric estimation of the ITF is as follows.(14)G^ωk,t=∑l=−Nb+NbG^lωkejlωsyst.For uncertainty analysis of ITF estimation and how to determine the optimal model order, see reference [29].


## 2.1. Periodic Time-Varying System

PTV systems are a class of common time-varying systems in practical engineering, such as electric motors, fans, shafts, and helicopter blades. The dynamic properties of these systems often contain a PTV coefficient.Since the number of teeth involved in meshing simultaneously in the gear transmission system changes periodically with time, the gear’s meshing stiffness, meshing damping, and other parameters will change periodically with time [27]. Therefore, a gear system can be considered a PTV system [28].Under the condition of uniform speed, the modeling and identification of the time-varying FRF of the gear system are carried out using the PTV system identification theory [29]with the torque generalized force and DTE as the input and output data of the system, respectively.


## 2.2. Input-Output Model of a Periodic Time-Varying System

The linear PTV systemG can be well described by a parallel structure consisting of infinitely many weighted linear time-invariant systems, as shown in Figure 2(a). In the figure, the triangle blocks denote a time-domain multiplication. Moreover, the (triangles) gains are written in a shorthand notation as lωsys instead of ejlωsyst. In practice, we use the branch structure composed of finite (such as Nb)weighted linear time-invariant systems to better approximate system G(1)Gω,t=∑l=−Nb+NbGlωejlωsyst,where Gω,t is called the instantaneous transfer function (ITF), Glω is called harmonic transfer function (HTF), and ωsys is the time-varying frequency of the system, namely, the meshing frequency of the gear transmission system.Figure 2
Principle diagram of a linear periodic time-varying system. (a) System direct equivalent model. (b) System indirect equivalent model.
(a)(b)In this paper, the ITF is identified in the frequency domain. The following are the basic assumptions needed in the identification.Without changing the input-output relationship of the system, we can exchange the time-varying gain module and the linear time-invariant module in Figure2(a) to obtain the model as shown in Figure 2(b). This model can be regarded as a linear time-invariant system with multiple inputs and a single output, and the input-output relationship in the frequency domain of the model can be expressed as follows.(2)Y0k=∑l=−Nb+NbHlωkU0k−pl,where U0k, Y0k are the spectra of discrete Fourier transform of input and output vectors. Hlωk is the HTF corresponding to the model in Figure 2(b), p=ωsys/ω0, and ωk=kω0.


## 2.3. Local Polynomial Estimation of the ITF

For the abovementioned model, this paper uses the local polynomial method (LPM). The LPM is a nonparametric identification method commonly used for linear time-invariant systems with multiple inputs and multiple outputs. This method uses polynomials to model FRF and transient terms in the local interval and uses the least square method to estimate model parameters based on smooth frequency characteristics [30].Assuming that the noise model is additive, according to (2), we can express the input-output relationship of a linear PTV system with output noise as follows.(3)Yk=HωkUk+WkUk=U0k.Uk is the input vector,(4)Uk=Uk−pNb…Uk…Uk+pNbT.The FRF between inputU0k and output Y0k is as follows.(5)Hωk=HNbωk…H0ωk…H−Nbωk.In addition,Wk is an error term consisting of the output noise (NY) and the output stochastic nonlinear disturbance (YS).In a local frequency band, that is, at the frequency linesk+l for (l=−n,−n+1,…,0,…,n), the FRF Hωk can be locally approximated by a low-order polynomial.(6)Hωk+l=Hωk+∑r=1Rhrklr+Ol/NR+1,where R is the order of the polynomial.Neglecting the bias term in (6), at the frequency line k+l, (3) can be expressed as follows.(7)Yk+l=ΘkKk+l+Vk+l,where(8)Θk=Hωk,h1k,h2k,…,hRk,is a 1×R+12Nb+1 row vector of the unknown parameters, and the vector Kk+l is defined as follows:(9)Kk+l=1l...lR⊗UTkT,where ⊗ is the Kronecker product operator. Evaluating (7) at l=−n,−n+1,…,0,…,n and grouping the equations in a matrix equation finally yield as follows:(10)Yk=ΘkKk+Vk,where Yk, Kk, and Vk are, respectively, 1×2n+1, R+12Nb+1×2n+1 and 1×2n+1 matrices.Through the least square method, we yield as follows.(11)H^ωk=YkSk,(12)S=KHKKH−1I2Nb+1,0H,where frequency argument k in (12) is removed for notational simplicity and xH is the complex conjugate transpose of x.The estimation of the HTFG^lω can be obtained from equations (5), (11), (13)(13)G^lωk=H^lωk+pl,where p=ωsys/ω0.Finally, the nonparametric estimation of the ITF is as follows.(14)G^ωk,t=∑l=−Nb+NbG^lωkejlωsyst.For uncertainty analysis of ITF estimation and how to determine the optimal model order, see reference [29].


## 3. Method Verification and Analysis

In order to verify the effectiveness and accuracy of the data-driven modeling method for the gear system. This paper establishes a three-degree-of-freedom gear-lumped mass model, and the time-varying FRF is identified.


### 3.1. Three Degrees of Freedom Gear Lumped Mass Model

The gear system is an elastic structural system with parametric excitation and nonlinearity. If the gear transmits a large load, the meshing surface of the gear is always in a contact state, and the backlash between the teeth will not affect the system’s dynamic characteristics [31]. This paper ignores the gear backlash, and the classical three-degree-of-freedom gear-lumped mass model is established. The model is more suitable for working conditions under heavy loads, as shown in Figure 3.Figure 3
Dynamics model of the gear system.Dynamic differential equations of the gear system with three degrees of freedom shown in Figure3 can be derived from Newton’s law.(15)mg1yg1″+cg1yg1′+chtp′t+kg1yg1+khtpt=0,mg2yg2″+cg2yg2′−chtp′t+kg2yg2−khtpt=0,mc1x″+chtp′t+khtpt=Fm+FaTt,.where ygi and θgi are the translational displacement and torsional vibration displacement of the i-th gear (i=1,2), respectively, Ig1 and Ig2 are the rotational inertia of the gears, mg1 and mg2 are the translational masses of the gears, dg1 and dg2 are the base diameters, pt is DTE, cg1Δcg2 and cht are the supporting damping and gear tooth meshing damping of the driving and driven gears, respectively, Tg1 and Tg2 are driving torque and load torque, respectively, kg1,kg2 and kht are the supporting stiffness of the driving and driven gears and the time-varying meshing stiffness, respectively.(16)xt=dg12θg1t−dg22θg2t,mc1=1dg12/4Ig1+dg22/4Ig2,where mc1 is the equivalent mass of the gear pair.(17)Fm=2Tg1mdg1=2Tg2mdg2,FaTt=mc1Tg1at2Ig1,where the torque (Tg1) is divided into average component (Tg1m) and variable component(Tg1at). Fm is the gear meshing force caused by the average component of torque, FaTt is the gear meshing force caused by the variable component of torque, and the sum of Fm and FaTt is torque-induced generalized force.Gear dynamic transmission error is as follows.(18)pt=xt+yg1t−yg2t−et,where et is static transmission error (STE). Generally, it is assumed that the STE changes periodically with the meshing frequency, and its Fourier series expansion is as follows:(19)et=∑i=1neicosiωt+ϕi,where ei is the i-order harmonic amplitude of static error. ω is meshing frequency, ϕi is the initial phase of the i -order harmonic.The potential energy method [14] is used to calculate kht, which is the periodic function of angular displacement, as shown in Figure 4.(20)cht=2ζkht4Ig1Ig2dg12Ig2+dg22Ig1,where ζ is the damping ratio.Figure 4
Gear meshing stiffness.cht is generally calculated by meshing stiffness kht.The abovementioned dynamic (15) can be converted to a state equation by introducing state variables x1=yg1, x2=y˙g1, x3=yg1, x4=y˙g2, x5=pt, x6=p˙t(21)X˙=AtX+BU,Y=CX+DU..The transfer function between generalized force caused by load torque and DTE can be obtained by the Laplace transform of state equation.(22)Gs,t=CsI−At−1B,where(23)At=010000−kg1mg1−cg1mg100−−khtmg1−chtmg100010000−kg2mg2−cg2mg2khtmg2chtmg2000001−kg1mg1−cg1mg1kg2mg2cg2mg2−khtmc1−khtmg1−khtmg2−chtmc1−chtmg1−chtmg2,B=000001mc1TC=000010D=0,where T is a transpose operator.


### 3.2. Model Setting

Torque variation componentTgla(t) adopts random multiline excitation. The relevant parameters are as follows: a sampling frequency fs = 215 Hz, a sampling point N = 217, periods P = 4, phase R = 100, amplitude A = 50 N m. Time-varying period Tsys = 0.5s, and the torque average component Tglm = 200 N m. The specific model parameters of the gear system are shown in Table 1.Table 1
Basic parameters of the gear system.
Parameters/UnitsDriving gearsDriven gearsGear numbers2047Modules/mm3Pressure angle/(°)20Gear width/mm20Contact ratio1.65Elastic modulus (Pa)2.06e11Poisson’s ratio0.269Supporting stiffness(N·m−1)3.5e9Supporting damping (Ns·m−1)48.661.6


### 3.3. Model Solving

In this section, the fourth-order Runge–Kutta method is used to solve the lumped mass model of the gear system. The Runge–Kutta method is an essential implicit or explicit iterative method for solving nonlinear ordinary differential equations.LetZ1=yg1,Z2=yg2,Z3=pt, the matrix form of (15) is (24),(24)MZ¨+CZZ˙+KZZ=PZ,where(25)M=mg1000mg20−mc1mc1mc1,CZ=cg10cht0cg2−cht00cht,KZ=kg10kht0kg2−kht00kht,ΔPZ=00Fm−mc1×e¨t+FaTt.(24) is transformed to obtain the state equation.(26)Z¨=−1MP−CZ×Z˙−KZ×Z.The fourth-order Runge–Kutta method is used to numerically solve (26).The model input is shown in Figure5(a), whose value is the sum of Fm and FaTt. The obtained model response is shown in Figure 5(b), see (18) for details. Since the input and output are periodic signals, only the results of a single period are shown in the figure.Figure 5
Input and output of three degrees of freedom gear lumped mass model. (a) Torque generalized force (b) Dynamic transmission error.
(a)(b)


### 3.4. Identification Results

This section considers the data information in the0∼8000 Hz frequency band, and the time-varying FRF identified by (14) is shown in Figure 6. It can be found that the FRF changes periodically with time, and the natural frequency corresponding to the single-tooth contact area is lower than that corresponding to the double-tooth contact area. This is because when the gear rotates, the alternate meshing of single and double teeth changes the inherent characteristics of the gear system. With the change of time, the first-order natural frequency and the second-order natural frequency periodically change in the range of 2385∼3011 Hz, and 6938∼7165 Hz, respectively.Figure 6
Frequency response function.At the same time, the theoretical time-varying transfer function is calculated by (22), as shown in Figure 7. Compared with the identification results of Figure 6, the time-varying FRF obtained by identification is in good agreement with the theoretical transfer function. They proved that the method can accurately capture the time-varying characteristics of the gear system.Figure 7
Theoretical frequency response function.


## 3.1. Three Degrees of Freedom Gear Lumped Mass Model

The gear system is an elastic structural system with parametric excitation and nonlinearity. If the gear transmits a large load, the meshing surface of the gear is always in a contact state, and the backlash between the teeth will not affect the system’s dynamic characteristics [31]. This paper ignores the gear backlash, and the classical three-degree-of-freedom gear-lumped mass model is established. The model is more suitable for working conditions under heavy loads, as shown in Figure 3.Figure 3
Dynamics model of the gear system.Dynamic differential equations of the gear system with three degrees of freedom shown in Figure3 can be derived from Newton’s law.(15)mg1yg1″+cg1yg1′+chtp′t+kg1yg1+khtpt=0,mg2yg2″+cg2yg2′−chtp′t+kg2yg2−khtpt=0,mc1x″+chtp′t+khtpt=Fm+FaTt,.where ygi and θgi are the translational displacement and torsional vibration displacement of the i-th gear (i=1,2), respectively, Ig1 and Ig2 are the rotational inertia of the gears, mg1 and mg2 are the translational masses of the gears, dg1 and dg2 are the base diameters, pt is DTE, cg1Δcg2 and cht are the supporting damping and gear tooth meshing damping of the driving and driven gears, respectively, Tg1 and Tg2 are driving torque and load torque, respectively, kg1,kg2 and kht are the supporting stiffness of the driving and driven gears and the time-varying meshing stiffness, respectively.(16)xt=dg12θg1t−dg22θg2t,mc1=1dg12/4Ig1+dg22/4Ig2,where mc1 is the equivalent mass of the gear pair.(17)Fm=2Tg1mdg1=2Tg2mdg2,FaTt=mc1Tg1at2Ig1,where the torque (Tg1) is divided into average component (Tg1m) and variable component(Tg1at). Fm is the gear meshing force caused by the average component of torque, FaTt is the gear meshing force caused by the variable component of torque, and the sum of Fm and FaTt is torque-induced generalized force.Gear dynamic transmission error is as follows.(18)pt=xt+yg1t−yg2t−et,where et is static transmission error (STE). Generally, it is assumed that the STE changes periodically with the meshing frequency, and its Fourier series expansion is as follows:(19)et=∑i=1neicosiωt+ϕi,where ei is the i-order harmonic amplitude of static error. ω is meshing frequency, ϕi is the initial phase of the i -order harmonic.The potential energy method [14] is used to calculate kht, which is the periodic function of angular displacement, as shown in Figure 4.(20)cht=2ζkht4Ig1Ig2dg12Ig2+dg22Ig1,where ζ is the damping ratio.Figure 4
Gear meshing stiffness.cht is generally calculated by meshing stiffness kht.The abovementioned dynamic (15) can be converted to a state equation by introducing state variables x1=yg1, x2=y˙g1, x3=yg1, x4=y˙g2, x5=pt, x6=p˙t(21)X˙=AtX+BU,Y=CX+DU..The transfer function between generalized force caused by load torque and DTE can be obtained by the Laplace transform of state equation.(22)Gs,t=CsI−At−1B,where(23)At=010000−kg1mg1−cg1mg100−−khtmg1−chtmg100010000−kg2mg2−cg2mg2khtmg2chtmg2000001−kg1mg1−cg1mg1kg2mg2cg2mg2−khtmc1−khtmg1−khtmg2−chtmc1−chtmg1−chtmg2,B=000001mc1TC=000010D=0,where T is a transpose operator.


## 3.2. Model Setting

Torque variation componentTgla(t) adopts random multiline excitation. The relevant parameters are as follows: a sampling frequency fs = 215 Hz, a sampling point N = 217, periods P = 4, phase R = 100, amplitude A = 50 N m. Time-varying period Tsys = 0.5s, and the torque average component Tglm = 200 N m. The specific model parameters of the gear system are shown in Table 1.Table 1
Basic parameters of the gear system.
Parameters/UnitsDriving gearsDriven gearsGear numbers2047Modules/mm3Pressure angle/(°)20Gear width/mm20Contact ratio1.65Elastic modulus (Pa)2.06e11Poisson’s ratio0.269Supporting stiffness(N·m−1)3.5e9Supporting damping (Ns·m−1)48.661.6


## 3.3. Model Solving

In this section, the fourth-order Runge–Kutta method is used to solve the lumped mass model of the gear system. The Runge–Kutta method is an essential implicit or explicit iterative method for solving nonlinear ordinary differential equations.LetZ1=yg1,Z2=yg2,Z3=pt, the matrix form of (15) is (24),(24)MZ¨+CZZ˙+KZZ=PZ,where(25)M=mg1000mg20−mc1mc1mc1,CZ=cg10cht0cg2−cht00cht,KZ=kg10kht0kg2−kht00kht,ΔPZ=00Fm−mc1×e¨t+FaTt.(24) is transformed to obtain the state equation.(26)Z¨=−1MP−CZ×Z˙−KZ×Z.The fourth-order Runge–Kutta method is used to numerically solve (26).The model input is shown in Figure5(a), whose value is the sum of Fm and FaTt. The obtained model response is shown in Figure 5(b), see (18) for details. Since the input and output are periodic signals, only the results of a single period are shown in the figure.Figure 5
Input and output of three degrees of freedom gear lumped mass model. (a) Torque generalized force (b) Dynamic transmission error.
(a)(b)


## 3.4. Identification Results

This section considers the data information in the0∼8000 Hz frequency band, and the time-varying FRF identified by (14) is shown in Figure 6. It can be found that the FRF changes periodically with time, and the natural frequency corresponding to the single-tooth contact area is lower than that corresponding to the double-tooth contact area. This is because when the gear rotates, the alternate meshing of single and double teeth changes the inherent characteristics of the gear system. With the change of time, the first-order natural frequency and the second-order natural frequency periodically change in the range of 2385∼3011 Hz, and 6938∼7165 Hz, respectively.Figure 6
Frequency response function.At the same time, the theoretical time-varying transfer function is calculated by (22), as shown in Figure 7. Compared with the identification results of Figure 6, the time-varying FRF obtained by identification is in good agreement with the theoretical transfer function. They proved that the method can accurately capture the time-varying characteristics of the gear system.Figure 7
Theoretical frequency response function.


## 4. Identification of FRF of Transient Dynamics Finite Element Model of Gear System

### 4.1. Establishment of a Transient Dynamic Finite Element Analysis Platform for Gear Systems

ANSYS software is a widely used and powerful finite element analysis software. The ANSYS Parametric Design Language (APDL) provides users with the integrated functions of parametric modeling, loading, solving, and postprocessing, which is very suitable for modeling and analyzing gear and other shape specification parts. The involute tooth profile of the gear can ensure a constant transmission ratio and constant positive pressure direction between tooth profiles [32]. According to the involute formation principle (Figure 8), the involute equation is established.(27)xθ=rbcosθ+θsinθ,yθ=rbsinθ−θcosθ,where rb is the radius of the base circle, θ is the expansion angle of k points on the involute.Figure 8
Involute formation principle diagram.In Figure8, rb is the radius of the gear base circle, αk is the pressure angle of the involute at k point, ra is the radius of the tooth top circle, and rf is the radius of the tooth root circle. The calculation formula is as follows:(28)rb=12mzcosα,ra=12mz+ha∗m,rf=12mz−hc∗+c∗m,where m is gear modulus, z is tooth number, α is pressure angle, ha∗ is top height coefficient，and c∗ is top clearance coefficient.The accurate geometric model is established using APDL, and the model is considered the elastic supporting for simulation analysis. To improve the efficiency and accuracy of finite element analysis; in this paper, when meshing the gear model, first, the involute cylindrical spur gear is divided by the sweeping method, and finally, make the local mesh of the gear tooth surface denser and finer, as shown in Figure9.Figure 9
Finite element analysis model of gear pair.Finite element modeling of the transient dynamics of gear systems, with the following considerations.(1)
Material: driving and driven gears have the same material.(2)
Contact pair: surface-to-surface contact. The tooth surface on the meshing side of the driving gears is considered the contact surface, and the tooth surface on the meshing side of the driven gear is considered the target surface. The most widely used penalty function, Lagrange’s algorithm, is chosen for the calculation.(3)
Reference point: gear failure is primarily gear tooth failure, so this paper simplifies the driving and driven gear’s hubs. The inner ring of the wheel hub is set as a rigidization circle. And the part between the rigidization circle and the gear rotation center is treated as a rigid body in the analysis process. To simulate the rotation of the gear around the axis, the guiding nodes are established at the center of the driving and driven gears, respectively, and the point-surface contact pair is established between the guiding node and the rigidization surface. The point-surface contact algorithm between the guiding node and the rigidization surface is defined as the multi-point constraint (MPC) algorithm. The outer ring of the wheel hub is set to be bonded to the gear’s inner surface, and the gear’s inner ring drives the gear rotation.(4)
Analysis step: we set the analysis type of ANSYS to transient analysis and specify the solution method for transient analysis as the default full method, including large deformation effects and a specific time step.(5)
Supporting set-up: the gear is installed on the gear shaft. Both of them are elastic, so the influence on the inherent characteristics of the gear body cannot be ignored. Especially when the gear shaft is thin and flexible, the difference will be more significant. This paper simplifies the complexity of the model and only considers the elastic support in they-direction. Therefore, the spring element is used to simulate the supporting stiffness in the y-direction equivalently, and the steering nodes of the inner ring of the driving and driven gears are coupled to the center mass point.(6)
Boundary conditions and loads: the driving and driven gears retain degrees of freedom for rotation about thez-axis and degrees of freedom for translation in the y-direction, constraints are applied in other directions, torque is applied to the driving gears, and rotational speed is applied to the driven gears.The angular displacement data on the driving and driven gears are extracted in ANSYS postprocessing, and the support reaction moment on the gear is extracted.


### 4.2. Model-Related Parameter Settings

The finite element model is intuitive and accurate, which can better reflect the actual working state of the gear system. Based on the virtual experimental data of the finite element model of the transient dynamics of the gear transmission system, this section uses the established data-driven modeling method of the gear system to identify the FRF. The relevant parameters of the finite element model of the gear transmission system are shown in Table2.Table 2
Basic parameters of the gear system.
Parameters/UnitsDriving gearsDriven gearsGear numbers2047Modules (mm)3Pressure angle (°)20Gear width (mm)26Contact ratio1.65Density (kg/m3)1410Elastic modulus (Pa)2.6 × 109Poisson’s ratio0.39The torque load variation component uses multisine excitation, and load data are generated in MATLAB software. The parameters are set as follows: a sampling frequencyfs  = 213 Hz, sampling points N  = 213, periods P = 4, phase R  = 1, amplitude A  = 0.4 N m. Time-varying period Tsys  = 0.5 s, the torque average component Tg1m  = 2 N m. The supporting stiffness kb  = 6 × 106 N m−1 and the supporting damping cb  = 8 × 104 Ns·m−1 of the gear system model when elastic supporting is considered.In the finite element software, the load data are read by the load file method. The torque load is applied to the driving gears, and the constant speed is applied to the driven gears for simulation.In ANSYS postprocessing, the angular displacement data of the driving and driven gears are extracted. According to (18), the DTE is generated, and the counter torque on the driven gear is extracted. The model excitation is shown in Figure 10(a), and the DTE is shown in Figure 10(b).Figure 10
Input and output of finite element model. (a) Torque generalized force. (b) Dynamic transmission error.
(a)(b)


### 4.3. Identification Results

This section considers the data information on the 0∼3000 Hz frequency band. It takes the torque generalized force and the DTE as the input and output data of the system, respectively. The time-varying FRF identified is shown in Figure11.Figure 11
Frequency response function.As shown in Figure11, with the change of time, the FRF changes periodically. The natural frequency corresponding to the single-tooth contact area is lower than that of the double tooth. This law is due to the change in the inherent characteristics of the gear system caused by the alternating meshing of single and double teeth when the gear rotates. With the change of time, the first-order natural frequency varies periodically in the range of 1005∼1556 Hz. This method is extended to the virtual experimental data of the gear system transient dynamics simulation platform, which still has high accuracy and specific engineering practical significance.


## 4.1. Establishment of a Transient Dynamic Finite Element Analysis Platform for Gear Systems

ANSYS software is a widely used and powerful finite element analysis software. The ANSYS Parametric Design Language (APDL) provides users with the integrated functions of parametric modeling, loading, solving, and postprocessing, which is very suitable for modeling and analyzing gear and other shape specification parts. The involute tooth profile of the gear can ensure a constant transmission ratio and constant positive pressure direction between tooth profiles [32]. According to the involute formation principle (Figure 8), the involute equation is established.(27)xθ=rbcosθ+θsinθ,yθ=rbsinθ−θcosθ,where rb is the radius of the base circle, θ is the expansion angle of k points on the involute.Figure 8
Involute formation principle diagram.In Figure8, rb is the radius of the gear base circle, αk is the pressure angle of the involute at k point, ra is the radius of the tooth top circle, and rf is the radius of the tooth root circle. The calculation formula is as follows:(28)rb=12mzcosα,ra=12mz+ha∗m,rf=12mz−hc∗+c∗m,where m is gear modulus, z is tooth number, α is pressure angle, ha∗ is top height coefficient，and c∗ is top clearance coefficient.The accurate geometric model is established using APDL, and the model is considered the elastic supporting for simulation analysis. To improve the efficiency and accuracy of finite element analysis; in this paper, when meshing the gear model, first, the involute cylindrical spur gear is divided by the sweeping method, and finally, make the local mesh of the gear tooth surface denser and finer, as shown in Figure9.Figure 9
Finite element analysis model of gear pair.Finite element modeling of the transient dynamics of gear systems, with the following considerations.(1)
Material: driving and driven gears have the same material.(2)
Contact pair: surface-to-surface contact. The tooth surface on the meshing side of the driving gears is considered the contact surface, and the tooth surface on the meshing side of the driven gear is considered the target surface. The most widely used penalty function, Lagrange’s algorithm, is chosen for the calculation.(3)
Reference point: gear failure is primarily gear tooth failure, so this paper simplifies the driving and driven gear’s hubs. The inner ring of the wheel hub is set as a rigidization circle. And the part between the rigidization circle and the gear rotation center is treated as a rigid body in the analysis process. To simulate the rotation of the gear around the axis, the guiding nodes are established at the center of the driving and driven gears, respectively, and the point-surface contact pair is established between the guiding node and the rigidization surface. The point-surface contact algorithm between the guiding node and the rigidization surface is defined as the multi-point constraint (MPC) algorithm. The outer ring of the wheel hub is set to be bonded to the gear’s inner surface, and the gear’s inner ring drives the gear rotation.(4)
Analysis step: we set the analysis type of ANSYS to transient analysis and specify the solution method for transient analysis as the default full method, including large deformation effects and a specific time step.(5)
Supporting set-up: the gear is installed on the gear shaft. Both of them are elastic, so the influence on the inherent characteristics of the gear body cannot be ignored. Especially when the gear shaft is thin and flexible, the difference will be more significant. This paper simplifies the complexity of the model and only considers the elastic support in they-direction. Therefore, the spring element is used to simulate the supporting stiffness in the y-direction equivalently, and the steering nodes of the inner ring of the driving and driven gears are coupled to the center mass point.(6)
Boundary conditions and loads: the driving and driven gears retain degrees of freedom for rotation about thez-axis and degrees of freedom for translation in the y-direction, constraints are applied in other directions, torque is applied to the driving gears, and rotational speed is applied to the driven gears.The angular displacement data on the driving and driven gears are extracted in ANSYS postprocessing, and the support reaction moment on the gear is extracted.


## 4.2. Model-Related Parameter Settings

The finite element model is intuitive and accurate, which can better reflect the actual working state of the gear system. Based on the virtual experimental data of the finite element model of the transient dynamics of the gear transmission system, this section uses the established data-driven modeling method of the gear system to identify the FRF. The relevant parameters of the finite element model of the gear transmission system are shown in Table2.Table 2
Basic parameters of the gear system.
Parameters/UnitsDriving gearsDriven gearsGear numbers2047Modules (mm)3Pressure angle (°)20Gear width (mm)26Contact ratio1.65Density (kg/m3)1410Elastic modulus (Pa)2.6 × 109Poisson’s ratio0.39The torque load variation component uses multisine excitation, and load data are generated in MATLAB software. The parameters are set as follows: a sampling frequencyfs  = 213 Hz, sampling points N  = 213, periods P = 4, phase R  = 1, amplitude A  = 0.4 N m. Time-varying period Tsys  = 0.5 s, the torque average component Tg1m  = 2 N m. The supporting stiffness kb  = 6 × 106 N m−1 and the supporting damping cb  = 8 × 104 Ns·m−1 of the gear system model when elastic supporting is considered.In the finite element software, the load data are read by the load file method. The torque load is applied to the driving gears, and the constant speed is applied to the driven gears for simulation.In ANSYS postprocessing, the angular displacement data of the driving and driven gears are extracted. According to (18), the DTE is generated, and the counter torque on the driven gear is extracted. The model excitation is shown in Figure 10(a), and the DTE is shown in Figure 10(b).Figure 10
Input and output of finite element model. (a) Torque generalized force. (b) Dynamic transmission error.
(a)(b)


## 4.3. Identification Results

This section considers the data information on the 0∼3000 Hz frequency band. It takes the torque generalized force and the DTE as the input and output data of the system, respectively. The time-varying FRF identified is shown in Figure11.Figure 11
Frequency response function.As shown in Figure11, with the change of time, the FRF changes periodically. The natural frequency corresponding to the single-tooth contact area is lower than that of the double tooth. This law is due to the change in the inherent characteristics of the gear system caused by the alternating meshing of single and double teeth when the gear rotates. With the change of time, the first-order natural frequency varies periodically in the range of 1005∼1556 Hz. This method is extended to the virtual experimental data of the gear system transient dynamics simulation platform, which still has high accuracy and specific engineering practical significance.


## 5. Conclusion

(1)
Based on the identification theory of PTV systems, this paper established the identification method of the time-varying FRF of the gear system. This method can accurately capture the change law of inherent characteristics with time, caused by the alternate meshing of single and double teeth during the rotation of the gear system.(2)
Furthermore, a virtual experimental platform of the gear system is established based on the finite element theory, and the data-driven modeling method is applied to the virtual platform. The identification results are also good, indicating that the method has specific practical significance in engineering.In addition, the method combined with speed measurement technology can be further extended to the frequency response function identification of gear transmission systems under variable speed conditions. It lays a foundation for fault diagnosis, parameter inversion, and vibration and noise reduction of the gear transmission system.



---

*Source: 1881361-2022-09-20.xml*