Dynamic Modeling of Gear System Based on 3D Finite Element Model and Its Application in Spalling Fault Analysis

A reduced-order dynamic model, based on three-dimensional (3D) finite element model (FEM) and component modal synthesis technique (CMS), was presented for simulating the dynamic behavior of the spur gear system. -e gear shaft and gear body were established via 3D elements to simulate bending and torsion of the gear system. -e CMS technique was used to generate a reduced-order model of a spur gear system. A pair of mating teeth was assimilated to two different foundations (one for the pinion tooth and the other one for the gear tooth) linked in series by some independent springs, which was used to simulate the contact stiffness. -e validity of the proposed model was verified by static analysis, dynamic analysis, and experimental analysis. -e results show that the proposed model is an effective model. In addition, the proposed model has also been applied to analyze spur gear spalling faults. -e results show that the dynamic response of the gear system is periodic vibration shock response due to the alternate meshing of single and double teeth. When the spalling fault occurs, some shock responses with significantly enhanced amplitude will be generated as the result of contact loss.


Introduction
Structural dynamical monitoring and fault diagnosis are of great importance to diagnose faults/failures of structures [1]. In order to facilitate the development of diagnostic and prognostic techniques for gear faults, many researchers worked on the gear dynamic modeling to ascertain the effect of gear defects on the vibration response. e first class of modeling method is the lumped parameter model. Ma et al. used a four degrees of freedom (DOF) gear dynamic model to study the dynamic and vibration characteristics of the gear system with local defects [2]. Chen et al. used a 6-DOF gear dynamic model to study the vibration of spur gear with tooth root crack propagating along tooth width and crack depth [3]. Kang et al. investigated the dynamic behaviors of the gearrotor system with viscoelastic supports under residual shaft bow effect [4]. Bartelmus et al. dealt with mathematical modeling and computer simulation as a tool for aiding gearboxes diagnostic inference [5]. Xiong et al. used single DOF model to investigate the effect of backlash on the nonlinear dynamic characteristics [6]. Wang et al. established a 3-DOF torsional vibration model and investigated the effects of the pinion speed and stiffness on the dynamic behavior of the gear transmission system [7]. Yang et al. investigated the effects of tooth breakage size and rotational speed on the vibration response of the planetary gearbox [8]. A simplified nonlinear lumped parameter model with 6-DOF established by Mohamed was developed to simulate the vibration response of faulty external spur gears [9]. Chen et al. developed a lumped parameter model to study the dynamic response of the gear system with different crack sizes [10]. e second approach is the hybrid beam element/lumped models. In this model, the beam elements are used to model the shafts and a pair of rigid disks is connected by a spring damper with time-varying mesh stiffness for modeling a pair of mating teeth. Wan et al. used this method to investigate the effect of tooth crack on the vibration response [11]. Han et al. analyzed the effects of position and type of the rotor crack (transverse or slant crack) on the gear system dynamic behaviors by this kind of dynamic model [12,13]. Ma et al. analyzed the fault feature of cracked gear [14,15]. Yu et al. used this model to study the dynamic coupling behavior of the transverse and rotational motions of gears subjected to gear eccentricities [16]. e first two models are relatively mature and have been used to analyze typical faults such as crack faults [17,18]. In addition, some researchers discussed the contact characteristic and gear geometry with existence of radial and angular eccentric errors [19]. Another model is the 3D finite element model. In recent years, the finite element modeling of some microstructure small cracks has made great progress [20][21][22]. However, it needs to be further studied in the equipment level system modeling, such as the gear system. Because the FEM is much computationally expensive compared with the first two models and mostly limited to quasi-static approaches [23,24]. To address this problem, a reduced-order gear dynamic model based on 3D FEM and CMS technique is presented and used to solve fast solution problems of 3D FEM. In addition, this model has also been applied to analyze spur gear spalling faults.
is paper is organized as follows. In Section 2, the 3D FEM of a spur gear system was developed. e Craig-Bampton component modal synthesis technique was used to generate a reduced-order model of the gear system. In Section 3, the tooth contact was modeled by using the time-varying elastic foundations (Pasternak foundations). In Section 4, the driving and driven gears are connected by the elastic foundation and the dynamic equations of the gear system are deduced. In Section 5, the equations of motion were solved by Newmark-β integration and the vibration characteristic of a gear system with/without spalling fault were simulated and discussed. Conclusions are given in Section 6.

The Gear Dynamic Model Based on
Craig-Bampton Component Modal Synthesis Technique Figure 1 shows the process of establishing the 3D FEM models of a spur gear system. Firstly, the FEM of the driving and the driven shaft is obtained by software HYPERMESH and ANSYS. e reduced-order model can be generated by Craig-Bampton component modal synthesis technique. en, the driving and the driven gear are connected by an elastic foundation with time-varying stiffness [23,25]. Finally, the dynamic equations are derived and solved by Newmark-β integration. e details of every step mentioned above are addressed in the following sections.

e Reduced-Order Model of Driving or Driven Gear.
A 3D FEM of driving shaft of a gear system is shown in Figure 2. ere are more than 3,000 nodes and 9,000 degrees of freedom in this model. It will be very time-consuming to solve the model directly. us, a reduced-order model is demanded for dynamic analysis. Component mode synthesis technique is a kind of efficient model reduction method. e Craig-Bampton (CB) component mode synthesis method is utilized in this paper.
As shown in Figure 2, according to the structure characteristic of the gear system, the 3D FEM of the driving shaft (including pinion and shaft) is divided into two substructures. ey are named as the shaft substructure and the gear substructure, respectively. e interface of the two substructures is the axle hole of the gear. e main applications of CMS on this model are briefly summarized as follows.
Assume that the physical coordinate of the driving shaft can be expressed as follows: where the subscripts "i"and "j" represent the internal node DOF and interface node DOF, respectively. "p" means stand for pinion, "k" denotes the kept DOFs used to connect the elastic foundation used to simulate tooth contact, andn the superscripts "s" and "g" represent the shaft substructure and the gear substructure. e mass and stiffness matrix of the shaft substructure was obtained by finite element analysis software and arranged in order: where M ii , M kk , and M jj are the mass matrix corresponding to internal node, kept node, and interface node, and other mass matrices are coupling matrices between different nodes. e meaning of the subscript of the stiffness matrix is the same as that of the mass matrix. By fixing all the interface coordinates, the branch characteristic equation of the fixed interface is obtained: where Φ s is normal modes and ω is Natural frequency. Additional constraints are introduced into all the interface DOF of the substructure, and then these interface DOF generate unit displacement in turn. e constraint modes can be obtained by solving the static balance problem: where R kk , R kj , R jk , and R jj are the reaction force when the unit displacement is applied on the interface. ψ s ik and ψ s ij are the residual constraint modes of the elastic foundation node and the fixed interface, respectively. e constraint modal set of substructure can be obtained by solving the following formula: e DOFs' transformation equations between physical coordinates and modal coordinates of the shaft can be written as follows: 2 Mathematical Problems in Engineering where φ s ir denotes a truncated set of normal modes Φ s and p s r , p s k , and p s j are the model coordinates corresponding to the mode mentioned above. en, the mass matrix, stiffness matrix, and force vector of the shaft substructure in the model coordinate can be expressed as follows: where M s p , K s p , and f s p denote the mass matrix, stiffness matrix, and force vector in the physical coordinate.
Similarly, the mass matrix M g pm , stiffness matrix K g pm , and force vector f g pm of the gear substructure in the model coordinate can be obtained by the same way. e removed fix-interface coordinates are synthesized by using interface displacement compatibility conditions: where S is the transfer matrix used to remove the nonindependent coordinate. e mass matrix, stiffness matrix, and force vector of the driving shaft in the independent model coordinate can be expressed as follows:  e same method is used to the driven shaft (including gear and shaft). e mass matrix M g , stiffness matrix K g , and force vector f g of the driven shaft can be obtained by the same way.

Mesh Stiffness Interface Modeling
In Section 2, the mass matrices, stiffness matrices, and force vectors of the driving and the driven shaft in the independent model coordinate were obtained. e relationship between them needs to be established by a pair of mating teeth. As shown in Figure 3, Bettaieb et al. [23,25] used the Pasternak foundation to simulate the tooth contacts. In their model, a pair of mating teeth is assimilated to two different foundations (one for the pinion tooth and the other one for the gear tooth). e two foundations are linked in series by some independent springs which was used to simulate the contact stiffness (see in Figure 3). In this paper, assuming that the two foundations are located in the plane of action, the stiffness of Hertzian contact of two meshing teeth can be given by [26] where E, L, and v represent Young's modulus, tooth width and Poisson's ratio, respectively. As shown in Figure 4, Pasternak's foundations are made of superposition of bending and shearing elements lying on independent springs. e differential equations of Pasternak's foundation can be written as follows [25]: where d(η) is the displacement at appoint of coordinate η; P(η) is the applied force; and D(η), G(η), and k(η) are the bending, shearing, and direct stiffness coefficients, respectively. As shown in Figure 5, assuming that the tooth face has n contact lines, every contact line is represented by a foundation.
en, the instantaneous strain energy of the ith foundation discretized into N cells is where k i,j , G i,j , and D i,j are the direct, shearing, and bending stiffness coefficients per unit of the length at the jth cell of the ith foundation; d i,j is the displacement of the foundation; δ i,j is the displacement of the gear body, which is corresponding to the kept DOFs p g k introduced in Section 2.1; Δη is the length of a cell; and k i,j can be deduced by the potential method [11]. e determination of G i,j and D i,j can be consulted in reference [25]. e stiffness of the ith foundation can be obtained by using finite difference and Lagrange's equation. Figure 6 is the integrated model of a gear system. A parametrically excited dynamic equation can be obtained by assembling all the mass matrices, the stiffness matrices, and the force vectors:

Gear Dynamic Modeling
where d is the displacement of elastic foundation; M qq K qq , and q are the mass matrix, stiffness matrix, and displacement  of the 3D FEM of the gear system, respectively; and K dd is the stiffness of elastic foundation and K dq is the coupling stiffness matrix of elastic foundation and 3D FEM of the gear system.
On the one hand, the mass matrix of equation (13) is a singular matrix. On the other hand, in order to reduce the number of DOF. e static condensation is needed for equation (13). According to equation (9), the foundation displacement vector d can be written as follows: en, equation (9) can be rewritten as follows: where K t,q � K qq − K qd K −1 dd K dq and f t,q � f q − K qd K −1 dd f dq . e damping model is the Rayleigh damping, C t,q � βK t,q , where β is a damping coefficient and β � 1 × 10 − 7 in this paper.

e Model of Spalling.
e references mentioned in the introduction mainly investigate the effect of gear tooth crack on vibration response. However, the spalling fault of gear is also a common failure type in engineering applications. e spalling fault usually happens near the pitch line of the gear. As shown in Figure 7(a), the spalling defect in this paper is modeled near the pitch line as a rectangular spalling. Figure 7(b) is the elastic foundation model corresponding to Figure 7(a). Assuming that the driving gear does not contact with the driven gear in the spalling area, the contact stiffness is removed in the spalling area. So, the contact stiffness in the spalling area can be obtained as follows: where ΔL is the reduction of the length of the contact line and other parameters have the same definition with equation (10).

Simulation and Experiment Results
e model validation is a very critical step. erefore, the validity of the model was verified by static and dynamic analysis. Figure 8, the node forces were applied on the whole tooth face and the three middle nodes in the pitch circle, respectively. Constraints were imposed on the nodes shown in Figure 8(a). e same forces and constraints are applied to the proposed model. Figures 9(a) and 9(b) are the node displacement across the face width when the force was applied on the whole tooth face and on the three middle nodes, respectively. Compared with FEM, the maximum relative errors of node displacement of the two cases are 3.93% and 4.72%, respectively. It indicates that the deflection curves from the proposed method are in well agreement with those results obtained by 3D FEM.

Dynamic Behavior Analysis.
In dynamic behavior analysis, two comparisons are used to verify the effectiveness of the proposed model. e first comparison is the comparison between the proposed method and experiments. e second comparison is the comparison between the proposed method and the hybrid beam element/lumped models.

e Comparison between the Proposed Method and
Experiments. Simulation and experimental parameters are listed in Table 1. A photograph of a spur gear test rig is shown in Figure 10. Simulation results of dynamic responses are obtained by Newmark-β numerical integration. e comparison results between experimental signals and the simulation signals are shown in Figures 11 and 12. Figure 11(a) is the time-domain waveform of the simulation signal. It shows that the system will produce impact response, when the single and double teeth meshing, alternately. Each set of impact is composed of two alternating periodic shock responses with different amplitudes. e frequency of every two strong shocks or every two weak Mathematical Problems in Engineering shocks is 394 Hz (0.0025 s), which is equal to the meshing frequency of the gear system. In the experimental signal, this phenomenon is not very clear due to the background noise. e second generation wavelet transforms is used to process the original signal, and the forth sub-band is shown in Figure 11(b). It shows some of the same phenomena as the simulation signal. Figure 12 shows the spectrum of simulation signals and experimental signals. It can be found that the vibration signals of the normal gear are mainly the meshing frequency and its doubling frequency.

e Comparison Results between Hybrid Beam Element/Lumped Models and the Proposed Method.
Based on the same parameters, the hybrid beam element/lumped models [11] and the proposed method were used to investigate the vibration characteristics, respectively. e comparison between them is shown in Figure 13. By comparison, it is found that the two models can show the impact response caused by the single and double teeth    is shows that the efficiency of the 3D dynamic model is solved. e validity of the model was also proved to some extent. Figure 14, assuming that the gear with the spalling fault is located in the driving gear, the spalling area is 2 mm long and 1 mm wide, and other parameters are shown in Table 1. Figure 15     Mathematical Problems in Engineering amplitude. e frequency of the shock response is equal to the rotation frequency of the driving shaft where the spalling gear is located. Figure 16 shows the spectrum of simulation signals and experimental signals of the gear with the spalling fault. By comparing Figures 12 and 16, it can be found that when the spalling fault occurs, the sidebands will appear in the spectrum. "Because of the short-period nature of the impact, the corresponding modulation sidebands will spread over a wide frequency range, which can produce high-order sidebands with low amplitudes. In the frequency domain, one can clearly observe sidebands with significant amplitudes starting in the resonance region (1446 Hz in Figure 16(a) and 1330 Hz in Figure 16(b)) and expand to higher frequencies [27]." ese characteristics were successfully captured by the proposed dynamic model.

Conclusions
Dynamic modeling can improve the understanding of the diagnostic information that is buried in the vibration signal. A reduced-order dynamic model based on 3D FEM and CMS is proposed in this paper. In addition, the proposed model was used to study the effects of tooth spalling on dynamic vibration response. Based on the materials presented in this paper, conclusions can be summarized as below: (1) e static analysis, dynamic analysis, and experimental analysis show that the proposed model is an effective model. e proposed model can be used to investigate the vibration characteristic of the spur gear system. (2) e simulation results and experiment results show that the system will produce impact response, when the single and double teeth meshing, alternately. Each set of impact is composed of two alternating periodic shock responses with different amplitudes. e frequency of every two strong shocks or every two weak shocks is meshing frequency.
(3) e gear spalling fault will cause shock response with significantly enhanced amplitude. e frequency of the shock response is equal to the rotation frequency of the driving shaft where the spalling gear is located.
Future research should consider other influences on the vibration response, such as the friction, the inter-tooth backlash, the vibration of housing, and manufacturing errors in gears.

Mathematical Problems in Engineering
Data Availability e data used to support the findings of this study are available from the corresponding author upon request.

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