Estimation of a Planetary Gear Mesh Stiffness: An Approach Based on Minimising Error Function

— The mesh stiffness of gear teeth is one of the major sources of excitation in gear systems. Many analytical and finite element methods have been proposed in order to determine the mesh stiffness of gears especially parallel axis spur gears. Most of these methods are not trivial because they involve complicated analyses which incorporate parameters like gear tooth error, gear spalling sizes and shapes, nonlinear contact stiffness and sliding friction before mesh stiffness can be determined. In this work, a method is proposed to determine the sun-planet and ring-planet mesh stiffnesses of a planetary gear system. This approach involves fitting a relationship between the measured natural frequencies from an experimental modal test and natural frequencies predicted using an analytical model of a planetary gear. This method is relatively easier compared to the existing methods which involve complicated analyses. For this study, the average mesh stiffness estimated is 12.5 MN/m.


I. INTRODUCTION
The determination of the mesh stiffness of planetary gear systems is critical to prediction of their dynamic responses. The accuracy of gear model (which is useful in prediction of its dynamic behaviour) highly depends on the whether the mesh stiffness is accurate or not. The mesh stiffness is the mesh force per deflection at the contact between the meshing teeth [1]. It can vary when the number of teeth in contact changes or when the load fluctuates. In some mathematical models, the gear mesh stiffness is assumed to be constant and used in a linear time invariant (LTI) system. Ozguven [2] classified the mathematical models developed in gear dynamics into five. For one of the groups, the lateral vibrations of the gear mounted on the shaft are considered in two mutually perpendicular directions such that there is whirling of the shaft. The torsional vibration of the system is usually considered along with the lateral vibrations. In this case, the stiffness of the gear teeth in the model is considered as either constant or time varying mesh stiffness. Cooley and Parker [3] in their paper discussed four modelling approaches. The one with constant mesh stiffness considers LTI equations. The dynamic excitation resulting from change in the mesh stiffness of the geared is sometimes approximated and represented by a quantity called static transmission error. Static transmission error is defined as the error in the motion transmitted due to elastic deformation of the teeth at the mesh.
A comparison of two approaches of calculating the mesh stiffness of spur gear was done by Cooley et al. [4]. In the first approach, the mesh stiffness was determined by dividing mesh force by mesh deflection. The second approach focuses on calculation of the local slope of the force-deflection curve about a nominal deflection. There is a significant difference between the two approaches. It shows that that the tooth profile modification does not change the stiffness amplitude for local slope calculations approach. There are differences in the mesh stiffness for gears with teeth profile modifications when using average slope approach.
In most literatures, tooth mesh stiffness is usually determined using computational models [4]. The fluctuation in mesh stiffness has been determined using finite element and contact mechanics [5]- [8]. Photoelectricity was used by Raghuwashi and Parrey [1] to measure the mesh stiffness of gear teeth in the presence of cracks. Keikbusch et al. [9] calculated combined torsional mesh stiffness of spur gears using two and three dimensional finite element (FE) models. The formula for the combined torsional mesh stiffness of spur gears was derived using the two dimensional FE model while the three dimensional FE models made it possible to carry out further studies on the helical gears.
Howard et al. [10] modelled the tooth stiffness using finite element analysis (FEA). Instead of modelling the bending stiffness, the static torsional mesh stiffness was considered and converted to a linear stiffness for use in their dynamic model. Chang et al. [11] proposed a model for determination of mesh stiffness of cylindrical gears using a combination of local contact analysis and finite element method (FEM). Sanchez et al. [12] presented a model for evaluation of spur gear mesh stiffness which includes parameters like bending, shear, compressive and contact deflections. They evaluated the meshing stiffness of spur gear pairs considering both the global and local tooth deformations. The stiffness was evaluated at any point of the path of contact by an analytical equation. The load at any point of contact can be determined from the analytical equation. Omar et al. [13] presented a method for calculating the gear mesh stiffness in the presence of large propagating crack in the root of a tooth. An analytical model for calculating the time-varying mesh stiffness (TVMS) of a spur gear was established by Hui et al. [14]. They considered the effect of extended tooth contact, revised fillet-foundation stiffness under tooth engagement region, nonlinear contact stiffness and tooth spalling defect in the analytical model. A method was developed by Fatih et al. [15] to determine the mesh stiffness of spur gears with asymmetric teeth. The first stage of their investigation focuses on calculation of a single tooth mesh stiffness using finite element method (FEM). The second stage focuses on calculation of gear mesh stiffness using already determined tooth stiffness result. Zaigang and Yimin proposed a general analytical model for calculating the mesh stiffness of spur gears [16]. They incorporated in the model, the effect of tooth errors which is as a result of tooth profile modifications, manufacturing errors or assembly errors. Ankur et al proposed an analytical method to calculate the TVMS of spur gear pair for different shapes, sizes and locations of spall considering sliding friction [17]. The result shows that these parameters are very important in calculating the TVMS of gears.
In this study, the mesh stiffness of sun-planet and ringplanet were determined using an error fitting method. Fitting method in this study involves minimisation of error in the squared difference between the predicted and measured natural frequencies as a function of mesh stiffness. This method which is relatively easier than the existing ones has not been used before to determine the mesh stiffness of a planetary gear. This was done by fitting a relationship between the natural frequencies measured from the experiment and those predicted by the analytical model. The model is similar to that of Lin and Parker [18].

II. THEORETICAL MODEL
The extended dynamic model of the planetary gear for this study is shown in Fig. 1. A fixed frame of reference was used to derive the equations of motion. The model comprises a carrier, ring, sun and two planet gears. The carrier, ring, sun and planet gears coordinates xh, yh, h=c,r,s and ζn, ηn, n =1, ….., N are with respect to the frame i,j and k, index n refers to the nth planet and N is the number of planet gears. The rotational coordinates are uj = rjθj. The stiffnesses kcx, krx, ksx and kcy, kry, ksy represent the bearing stiffnesses of the carrier, ring gear and sun gear in x and y coordinates, krn and ksn are the ring-planet and sun-planet mesh stiffnesses respectively. The bearing stiffness of the planet gear is denoted by kp. The torsional stiffness of the carrier, ring, sun, and planets are denoted by kcu, kru, ksu and kpu respectively. The ring gear is assumed to be stationary; therefore, its torsional stiffness is very high. The torsional stiffness of the carrier, sun and planet gears are set to zero.
The following assumptions were made in formulating the equations of motion • The main components (i.e., the carrier, ring, sun, and planets) are all assumed to be rigid bodies while the connections (i.e., meshing) between sun-planet and ringplanet are represented by linear springs acting along the pressure line which is tangential to the base circle. Timevarying mesh stiffness due to the fluctuation of the number of tooth pairs in contact is ignored.
• The bearings are assumed to be represented by two perpendicular linear springs for each gear component. They are also assumed to be isotropic.
• The frictional forces between the meshing teeth are neglected because the system is assumed to be well lubricated.
• The gear is free of manufacturing, profile, and other form of errors; therefore, there is no transmission error for the case of free vibration analysis.
• The two planet gears are identical and the angle between them is 180 o .
• The natural frequencies are determined for a linear timeinvariant case. The sun-planet mesh which was modelled by a linear spring is shown in Fig. 2. The subscript n refers to the nth planet. It shows the sun-planet mesh model which shows the pressure angle, (which is the angle between the pressure line AC and the common normal EF to the two gear teeth at the point of contact).
Then the mesh deflection at the sun-planet mesh is: Similarly, the mesh deflection at the ring-planet mesh is: Using Newton second law of motion, the equations of motion for the sun gear are written as: The general form of the equation in matrix form can be written as: where, , , ℜ × ; = , = , = , = are the mass, damping, bearing and contact stiffnesses matrices, respectively. The mass and stiffness matrices of the carrier, ring, sun and planet gears can be seen in the Appendix. Although two planet gears are shown in Fig.  1, the equation of motion can be used for a planetary gear with more than two planet gears. The parameters of the planetary gear system are shown in Table I. The predicted natural frequencies were determined using eigenvalue problem written as: where and ∅ are the predicted natural frequencies and vibration modes respectively. This was done for a case where the ring gear was excluded from the model so that the natural frequencies can be determined when only the sun-planet mesh is included. Ring gear was subsequently included in the model and the natural frequencies of the system were determined. Some of the predicted natural frequencies for the cases are shown in Tables II and III. III. EXPERIMENTAL MODAL TESTING In this section, experiments were carried out to measure the natural frequencies of the planetary gear system. Two cases were considered. For the first case, the ring gear was excluded from the test rig during the experiment. Therefore, only the carrier, sun and planet gears were present. For the second case, ring gear was included in the experiment. An instrumented hammer (model number PCB 086CO3) was used to excite the components. The sensitivity of the hammer was 2.25 mV/N. The response was measured by accelerometers mounted on either the carrier, ring, sun and planet gears in x and y directions. The sensitivity of the accelerometer was 1.00 mV/ms -2 . The sampling rate was 2.56×10 kHz. Signal analyser (Fig. 3) was used to obtain the frequency response functions (FRFs) and the Nyquist circle. Only the Nyquist circle was shown in Fig. 4 and the natural frequencies were recorded. The measured natural frequencies for case 1 are shown in Table II. They were used in the model to determine the sunplanet mesh stiffness.
Likewise, the measured natural frequencies for the case two are shown in Table III, they were subsequently used to determine the ring-planet mesh stiffness. For instance, the Nyquist circle when a natural frequency of 515 Hz was measured is shown in Fig. 4.

IV. ESTIMATION OF MESH STIFFNESS BY FITTING METHOD
This section focuses on estimation of mesh stiffnesses by minimising the error function. The contact stiffness, Km matrix (shown in the Appendix) was decomposed into two matrices such that the sun-planet mesh stiffness can be varied. A reasonable range of estimated values for the mesh stiffnesses were assumed between 10 6 and 10 9 N/m. The relationship in (10) was used to determine the mesh stiffness by choosing to minimise an error function ∈.
where, n is the chosen number of modes to be used in the estimation process, while ωie and ωip are the natural frequencies of ith mode from the experiment and the predictions respectively. Table II shows both the chosen measured and predicted natural frequencies. Each pair of measured and predicted natural frequencies was chosen based on the closeness in their values. Some of the measured natural frequencies cannot be identified by the analytical model because of the assumptions in formulating the model which can only capture the natural frequencies of the translational and rotational modes.  There are two different lines of the error function ∈ for sunplanet mesh stiffness shown in Fig. 5. The red line denotes the error function when four natural frequencies were chosen from Table II while the black is for the error function when three natural frequencies were chosen. It shows that the same value of mesh stiffness was estimated whether three or four natural frequencies were considered. The minimum value in the error function corresponds to an estimated mesh stiffness of 1.0 × 10 7 N/m and this was taken as the sun-planet mesh stiffness in the analytical model. The same procedure follows for estimating the sun-gear mesh stiffness was repeated when a ring gear was included on the test rig and the mathematical model in order to estimate ring-planet mesh stiffness. The estimated ring-planet mesh stiffness as shown in Fig. 6 is 1.5 × 10 7 N/m. The average of sun-planet and planet-ring mesh stiffnesses estimated is 1.25 × 10 7 N/m. This value was subsequently taken as the estimated mesh stiffness for both the sun-planet and planetring contact. Since the same numbers of pairs of teeth are in contact at any point in time for the sun-planet and planet-ring mesh, both mesh stiffnesses should be the same provided the gear material are the same.

V. CONCLUSION
The mesh stiffness of a planetary gear has been estimated using an error fitting method. This is a new approach by which sun-planet and ring-planet mesh stiffnesses of a planetary gear system can be estimated separately. The average mesh stiffness is then calculated and taken as the time-invariant mesh stiffness. The mesh stiffnesses for the sun-planet and ring-planet are the same because the number of pair of teeth in contact are the same for both cases (i.e. they have the same contact ratio). In addition, the material selected for the sun, ring and planet gears are the same. The method involves minimisation of error between the natural frequencies obtained from an analytical model and experimental modal test. It is relatively easier compared to the existing complicated analytical models formulated by gear researchers to determine the mesh stiffness. For example, there is no need for any complicated teeth contact analysis as seen in literatures. This method can also be used to estimate the mesh stiffnesses of gears in contact for simple and compound parallel axis gears. The value of the estimated mesh stiffness obtained can be used in a mathematical model of a gear system in order to predict the dynamic behaviour. The Km matrix can be decomposed into three matrices such that the sun-planet and planet-ring mesh stiffness can be varied. The second matrix term was set to zero for sun-planet mesh stiffness estimation.

ACKNOWLEDGMENT
The author would like to appreciate Tertiary Education Trust Fund (TETFUND), Nigeria for sponsoring this research.