Dynamic Modeling of Gears: An Innovative Hybrid FEM–Analytical Approach

Gearboxes are mechanical systems characterized by special failure modes such as tooth root bending and surface fatigue (pitting). While the first is very dangerous and leads to cata- strophic failures of the whole system, the latter is less critical and related to a progressive degradation of the Noise, Vibration and Harshness (NVH) behavior that affects the perfor- mances of the system, before evolving into more severe consequences.

Moreover, it should be mentioned that in specific applications like wind turbines, gear- boxes are often the weakest components (major failures estimated generally every 5 years [1]). The capability to model the real-time behavior of geared transmissions under operation is very useful to plan maintenance interventions reducing the costs and preventing further damages.

Austin [2] studied, by means of numerical approaches, the behavior of a multi-stage trans- mission. Based on a test rig owned by NREL, he validates his models. Guo et al. [3] studied the load sharing in a 750-kW wind turbine gearbox. Kahraman et al. [4] evaluated the influence of the rim gear size on the deformation of the gearbox. The first results were promising, even though some discrepancies were observed. Experiments were also performed by Qu et al. [5]. They used fiber Bragg gratings and piezoelectric strain sensors mounted near the root of the tooth to measure the effect of pits on the gearbox behavior. Differences in the measures were observed between damaged and undamaged gear flanks. The knowledge of the behavior of the gearbox under different damage conditions might be useful to set up a sensor-based approach to understand the damage level which would help in planning maintenance interventions.

The aim of this preliminary work is to validate the hybrid finite elements (FE)–analytical approach on one side and to determine whether the presence and evolution of a surface damage (pitting) can be evaluated from the numerically predicted vibrational spectra on the other side.

Therefore, the preliminary activity consisted in the definition and validation of an approach that is able to model the effect of surface damages. At this stage, healthy gearboxes were modeled and the results compared with measurements. The numerical analysis of the selected geometries (a back-to-back test rig and a planetary gearbox) already represents a hard engi- neering challenge.

Once the models were validated with experimental data, the subsequent step was the sim- ulations of the system in the presence of different levels of damage to understand their effect on the NVH.

Two completely different geometries were used for validating the approach: a back-to-back test rig and an industrial planetary gearbox.

A back-to-back test rig was the first geometry used (Fig. 1) for validating the numerical approach. This test rig, available in the labs of Politecnico di Milano, is made by two parallel axes gearboxes connected by mean of two shafts forming a closed mechanical loop. The gearboxes have the same gear ratio 17/18, but a different number of teeth (34/36 and 17/18, respectively). The 34/36-teeth gearbox (called service) has helical gears, while the second one (called testing) has spur gears. The torque is applied by means of a relative rotation of the two semi-shafts (Fig. 1) that connect the two gearboxes (by means of a hydraulic actuator). An e-Motor supplies only the power losses (mainly due to friction) and not the entire power that flows in the mechanical loop.

Each gear of the service gearbox is supported by means of four tapered roller bearings (Fig. 2), while the gears belonging to the test gearbox are supported by two roller bearings.

The second gearbox modeled is a single-stage planetary gearbox having a reduction ratio i = 3. It is a low backlash precision planetary gearbox suitable for mechatronic applications. The sun gear is force fitted on the input shaft (Fig. 3) that is supported by means of a ball bearing. The planet carrier (output shaft) is supported, in turn, by means of two ball bearings. Three planets engage with sun gear and ring. The latter is directly broached in the ductile iron [6, 7] casted housing. The planets are supported by means of full-complement needle bearing [8, 9]. The complete data are reported in literature [10–13].

The modeling of the above-mentioned test rigs was made taking advantage of an innovative hybrid FE–analytical approach. The simulations were performed with the software Transmis- sion3D by ANSOL. It has been specifically developed to perform simulations of entire transmissions in a reasonable amount of time. It allows to model gears, bearings (with all the rolling elements), shafts, housing and the other typical machine elements involved in a gear transmission.

While a traditional FE approach is used to describe the macroscopic deflection of the com- ponents, the contacts between bodies (i.e. gear teeth, rolling elements and races in bearings, etc.) are solved analytically using the Hertzian theory.

A traditional numerical approach requires, in order to properly simulate a contact, a mesh refinement to ensure the convergency. Considering the complexity of a geared system and the fact that the position of the contact as well as the number of engaging surfaces is variable during operation, with a traditional FEM, each gear flank must be discretized with a very fine mesh, leading to an enormous increase in the number of equations to be solved. On the con- trary, the present approach solves a relatively coarse grid (few equations) in addition to an analytical equation for each active contact.

In detail, the points at which the individual contacting surfaces belonging to the two bodies are closest to each other before the application of load are found [14]. Subsequently, the principal normal curvatures of the surfaces are determined. With this information, the size of the contact zone is estimated using the Hertzian theory and a grid (larger than the predicted Hertzian contact ellipse) is laid out around each principal contact point (Fig. 4). The grid is then projected on the body surfaces. The computed contact pressures are not very sensitive to the size of the grid.

The displacement u(rij ; r) of a field point r due to a load at the surface grid point rij can be expressed as

where $\mathbf{q q}$ is a point inside the solid body, sufficiently far from the surface (Fig. 5). The first term, $\left[\mathrm{u}\left(\mathbf{r}_{\mathrm{ij}} ; \mathbf{r}\right)\right.$ and $\left.\mathrm{u}\left(\mathbf{r}_{\mathrm{ij}} ; \mathbf{q}\right)\right]$, is evaluated using the surface integral approach. The second one, $\mathrm{u}\left(\mathbf{r}_{\mathrm{ij}} ; \mathbf{q}\right)$, is obtained from the FE model [14].

The term in square brackets represents the deflection of $\mathbf{r}$ with respect to the point $\mathbf{q}$. This relative component can be better estimated using the Boussinesq half-space solution [15] rather than using the FE results. The deformation of the body will, in fact, not significantly affect this term. On the other side, the remaining term $\mathrm{u}\left(\mathbf{r}_{\mathrm{ij}} ; \mathbf{q}\right)$ is not significantly affected by local stresses at the surface. Thus, if $\mathbf{q}$ is far enough from the surface, $\mathrm{u}\left(\mathbf{r}_{\mathrm{ij}} ; \mathbf{q}\right)$ is better estimated from the FE model of the solid body. The location of $\mathbf{q}$ is called the 'matching' point. In order to match the surface integral and FE solutions, a set of points $\Gamma$ can be used instead of a single point $\mathbf{q}$ (Fig. 5) [16]. Additional details are given in [17].

While this approach has proven to be significantly more computationally efficient with respect to traditional FEMs, further improvements of the models are possible: an even leaner model was obtained by substituting the bearings with equivalent radial springs having an equivalent stiffness obtained from preliminary steady-state simulations.

The planetary configuration used in this study was also modeled with a traditional FEM approach. A comparison made between the computational effort required by the traditional and the actual approach on the same hardware (115GFLOPS) showed that the latter is about 700 times faster if the bearings were modeled and 2000 times faster if equivalent springs were used. While a quasi-static simulation with the traditional approach requires about 30 h, the simulation time of a single time step is less than a minute with the hybrid approach.

The tests on the back-to-back test rig were performed at a rotational speed of 3000 rpm $\left(\mathrm{GMF}_{17 / 18}=850 \mathrm{~Hz}\right.$ and $\left.\mathrm{GMF}_{34 / 36}=1700 \mathrm{~Hz}\right)$ without the applied load. An accelerometer mounted on the test gearbox allows the measurements of vibrations during operation (Fig. 1, left).

The planetary architecture was tested on an electrical closed-loop test rig made by an electric motor, a generator and two inverters. The meshing of the teeth $(\mathrm{GMF}=1800 \mathrm{~Hz})$ was the main source of excitation. Vibrations were measured by means of two piezoelectric accelerometers (PCB M352C68) fixed on the housing of the gearbox by means of magnetic supports (Fig. 2, left).

Table 1 shows the comparison between the numerical and experimental data for the back-to- back test rig and for the planetary gearbox.

For the back-to-back architecture, the main peaks appear at around 850 and 1700 Hz (meshing frequencies of the two gear pairs). The other peaks above the GMF (2550, 3400, etc.) are harmonics of the meshing frequencies. The first peak around 600 Hz is an eigenfre- quency of the rig. The agreement between the predicted and measured frequencies is good, except for this first frequency where a discrepancy of about 15% was observed. This is prob- ably due to imperfect modeling of the inertia of the e-Motor.

Also, for the planetary gearbox, comparison of the excited frequencies showed the good- ness of the numerical model.

After the approach was validated, surface damages were introduced in the numerical models.

In order to understand the effect of a surface damage on the vibrational spectrum, additional

simulations were performed on both gearbox configurations, including damages differing in

severity and the position on the flank. In this stage, the damage was supposed to be located

on one single teeth. The presence of damages produces a modification of the contact pattern.

As an example, the contact pressures resulting on a tooth flank are shown in Fig. 6 both for a healthy and a damaged flank.

The upper contour in Fig. 6 corresponds to a non-damaged tooth. The upper stressed region is due to the gear contact, while the one at the bottom is related to the bending stress. The stresses are quite uniform, and the small differences along the axial direction are due to deflection of the planet carrier. The lower contour in Fig. 6, which corresponds to a damaged tooth, shows a central region which is not significantly loaded, while the remaining portion of the flank is overloaded.

Even if the state of stress (and strain) could represent a potential indicator of the damage, in this initial phase, only vibrations that are easier to measure have been considered as a trig- ger for the damage initiation.

Figures 7 and 8 show the results of the numerical transient simulations ofthe back-to-back test rig for two different levels of damage and the healthy condition in terms of accelerations (vibrations). Results are reported both in the time as well as in the frequency domain.

The spectra show that the presence of pitting does not modulate the meshing frequencies and their harmonics: nevertheless, the nucleation and progression of the damage promote the excitation of additional frequencies.

This can be an important feature in the development of a model-based SHM algorithm for the assessment of the health of gearboxes [18–21], but some additional activities are neces- sary, in order to better define the quantitative relationship between the severity of the damage and the effect produced and to define the corresponding thresholds.

The validation of the capability of the model to predict the excited frequencies represents the goal of this preliminary step.

A new hybrid FEM–analytical approach has been proposed for the simulation of complex gear sets. Initial simulations of the real configuration – i.e. including bearings – were per- formed in order to calculate the stiffness of the supports. Subsequently, a leaner model in which all the supports were substituted by means of equivalent radial springs was used to model a back-to-back test rig and an industrial planetary gearbox.

The results of the transient simulations in terms of vibrational spectra were compared with the experimental measurements for both configurations and they showed a good agreement.

For the back-to-back test rig, the effect of the presence of pits on the vibrational spectra was investigated using a numerical approach: the simulations pointed out that the defects promote the excitation of additional frequencies (mainly below the GMF). This phenomenon increases with the severity of the damage. Nevertheless, a quantitative relation has still to be determined. Concerning the planetary gear set, the results of the first simulations in healthy conditions are promising, but still to be refined, and for this reason, the subsequent simula- tions in damaged conditions have been delayed.

Further studies are ongoing in order to complete the activity, with the aim of fine-tuning the numerical models and to better understand the quantitative relation between damage and frequency spectra. The final objective is the extension set up an SHM strategy for gears.