A New Analytical Method for Modeling the Effect of Assembly Errors on a Motor-Gearbox System

The well-known gear tooth defects such as root cracks and flank spalls have been widely investigated in previous studies to model their effects on the time varying mesh stiffness (TVMS) and consequently the dynamic response of motor-gearbox systems. Nevertheless, the effect of assembly errors such as the center distance and the eccentricity has been less considered in past works. Determining the signature of these errors on the system response can help for their early detection and diagnostic to avoid overloading and failure of gears. An original geometric-based method combined with the potential energy method is proposed in this paper to accurately model the effect of these assembly errors on the TVMS of mating spur gear pairs. This is achieved by updating the line of action equation (LOA) at each meshing step using the actual coordinates of gear centers and employing a contact detection algorithm (CDA) to determine the actual contact points coordinates. An electrical model of a three-phase induction machine was then coupled with a dynamic model of a one-stage spur gear system to simulate the effect of assembly errors on the electromechanical response of the motor-gearbox system. The simulation results showed that the center distance error induces a reduction in the TVMS magnitude and the contact ratio, whereas the eccentricity error causes a double modulation of the TVMS magnitude and frequency. In addition, the results showed that assembly errors can be detected and diagnosed by analyzing the system vibration and the motor phase-current.


Introduction
Thanks to their robustness and good efficiency-cost ratio, induction motors and gearboxes are widely adopted in modern power transmission applications. Gearboxes offer a large variety of solutions to realize customized power transmission needs in a very reduced space by choosing the convenient gear types (spur, helical, bevel, etc.), gears configuration (parallel shaft, planetary, etc.) and the right combination of teeth and stage numbers depending on the desired transmission ratio. To ensure a smooth power transmission and a uniform load distribution, gears should be properly designed and assembled in their housing [1]. Otherwise, design and mounting failures lead to the apparition of transmission problems such as backlash and excessive noise. Two frequent problems are almost inevitable when assembling gear systems: the center distance and the eccentricity errors.
The center distance error occurs when the distance separating the two mating gear centers is not adjusted properly, whereas the eccentricity error appears when a shift exists between the geometric center of the gear and its rotation center. These anomalies should be avoided in gear systems because they affect the transmission efficiency of gearboxes and reduce their practical useful life by increasing the tooth surface contact stress and the tooth root bending stress [2]. Understanding the effect of the center distance and eccentricity errors on the system's observable quantities such as the vibration and the electrical response could help enormously in the early detection and reparation of them in order to avoid system failure [3,4].
Dynamic modeling has always been an efficient tool to predict the response of motorgearbox systems under varying operating conditions. A set of models were proposed in the literature to simulate gearboxes and induction machine responses; among them, numerical models based on the finite element method (FEM) and analytical models based on the lumped parameter method (LPM) were the most popular ones [5][6][7][8][9]. FEM-based models generally provide more precise results and a high fidelity to the original system. However, they are very sensitive to the mesh density and the types of finite elements selected and their computational cost varies exponentially with the mesh density. In the LPM-based models, the system components are replaced with their respective inertial or electrical properties in concentrated points, and the interactions between them are modeled with basic physics laws such as the Newton's second law or the Maxwell-Faraday equations of electromagnetic induction. Analytical methods make it possible to obtain very satisfying results of the system performance and offer a good precision-complexity ratio.
Another method used to model the effect of mechanical faults on the induction motor current signature can be found in [10,11]. In this method, the different characteristic frequencies such as gear mesh, gear faults, etc., are directly introduced as fluctuations around the average load torque with subjectively chosen amplitudes based on experimental observations. Then, this torque is used as input of the electrical model of the induction machine to evaluate its effect on the stator phase current. This method allows us to easily model the electromechanical system, but it fails in quantifying the system excitations, especially the fault oscillation components.
A key parameter in modeling gear dynamics with the LPM is the gear mesh stiffness. Stiffness of mechanical structures is an inherent parameter that depends mainly on the structure geometry and the loading conditions; the gear mesh stiffness reflects the teeth capacity to support the contact load during meshing process. Various methods were proposed to estimate this parameter experimentally [12], numerically [13] or analytically [14]. Experimental and numerical methods for the TVMS calculation are mainly based on the determination of the gear deflection when subjected to the applied load, whereas analytical methods use the different energies stored in the gear tooth body when supporting the load to derive the gear mesh stiffness components by the potential energy method proposed in [15].
In the last decade, countless works have been published in the literature to investigate the effect of local gear teeth faults such as cracks and spalls on the gear mesh stiffness [16][17][18][19]; however, only few studies were realized to model the center distance and eccentricity errors effect on this parameter [20]. In [21], the authors studied the effect of eccentricity errors on the dynamic behaviour of planetary gears, where the time varying mesh stiffness (TVMS) was calculated using a simplified analytical formula in terms of elastic properties. A kinematic model based on the potential energy method was proposed in [22] to determine the effect of time varying center distance on the TVMS curve for perfect gears and gears with local gear teeth faults (cracks and spalls). In [23], planetary gear train motions and the mesh stiffness were solved simultaneously considering the time varying center distance caused by the eccentricity errors. In [24], the effect of eccentricity error on the TVMS and the dynamic response of a two-stage spur gear system was determined using the potential energy method and the lumped parameter theory. The effect of eccentricity error on helical gears was calculated in [25].
Experimental and FE methods for the TVMS determination need specific tools and materials and are often time and cost consuming, which is why significant efforts have been made to develop efficient analytical methods of calculating the gear mesh stiffness. Traditional analytical methods were based on the angular expressions of the TVMS components to determine the contact points, but these expressions make the reproduction of gears varying geometries and defects a real challenge. In addition, these mathematical expressions need huge transformations in order to consider gear center displacement during the meshing process.
An original geometric method based on the potential energy method is proposed in this paper to accurately model the effect of gears center distance and the eccentricity errors in the TVMS of mating gear pairs. This method was first proposed and validated in [26] to determine the effect of teeth surface defects (spalls and pits) on the TVMS and is further improved in this study to consider gear center displacement and varying teeth shapes. This is achieved by updating the line of action (LOA) at each meshing step using the actual coordinates of gear centers and employing a contact detection algorithm to determine these contact point coordinates. An electromechanical model of a motor-gearbox system is then developed to simulate the system dynamical and electrical response by combining an electrical model of the three-phase induction machine with the dynamical model of a one-stage spur gear system. The gearbox model was improved to consider the presence of assembly errors in the system.
The proposed method for the TVMS calculation presents many advantages compared to traditional analytical methods: • It uses linear expressions instead of angular expressions for the evaluation of the TVMS components, which makes the incorporation of gear faults less complicated; • This method could be seen as a compromise between the FEM and the analytical method, because the discretization of the gear allows one to consider the actual geometry as in the FEM, and then the equations derived from the potential energy method are used for the calculation of the TVMS because they are time efficient compared to the FEM; • It calculates the gear mesh stiffness based on the actual contact points and allows us to visualize the meshing process and stiffness evolution; • All modifications made to the gear tooth geometry could be incorporated in the TVMS calculation process such as tooth tip corner, the tooth width variation and the tooth root geometry.
This paper is organized as follows: Section 2 describes the method proposed in this study for motor-gearbox modeling, contact points detection, geometric errors incorporation and the stiffness calculation process. In Section 3, the simulation results are presented and discussed for different stages of the center distance error and different configurations of the eccentricity error. Finally, we conclude the paper with some summarizing conclusions. Figure 1 represents the overall flowchart of the TVMS evaluation method proposed in this study. The gears used in this investigation are industrial KHK gears of references: KHK SS2-29/KHK SS2-36 (these two gears were chosen because we have them in our laboratory, so we can verify the extracted geometry by the real gears, but the developed method is applicable for any type of gear). In the initialization part of this method, the gear 3D models downloaded from the constructor website are incorporated in a computer-aided design software to extract the gear teeth geometric features (the tooth profile, the tooth width variation, etc.). These features are used then to reconstruct the gear matrices as explained in Section 2.1.1. In the calculation part or the dynamic part of this flowchart, these matrices are transformed at each simulation step to adjust the position and the angular rotation of their centers; then, the corresponding LOA equation is defined and swept to detect contact points at every iteration. Finally, these points are used to calculate the stiffness at this iteration using Equations (1)-(5) and the angular position of the pinion is incremented for the next iteration. The TVMS evaluation steps are described in more detail in the following subsections.

Gear Representation in the Proposed Method
Each gear is modeled with an N × Z matrix as illustrated in Figure 2, where N represents the sampling rate of the chosen gear tooth profile and Z represents the number of gear teeth. Every cell (n, z) of this matrix contains the coordinates of the nth point of the zth tooth profile P z,n as shown in Figure 2. We should note here that only the active side of the gear teeth is considered in this representation.
This representation offers a set of advantages compared to the traditional analytical methods, because it allows one to efficiently consider all types of geometrical modifications or errors of the gear teeth. These modifications include: form errors (profile form error, lead form error, etc.) and location errors (single and multiple pitch errors, runout errors, etc.). Global tooth profile errors/modifications are incorporated directly in the tooth profile vector T P , whereas local profile errors are assigned to the specific tooth index. Runout error could be considered by performing a periodic translation of teeth profile vectors T P .

Gear 2D model
Gear tooth discretization Gear matrix The above mentioned errors could be handled by the two-dimensional representation proposed in this study; other high dimensional modifications such as gear tooth crowning, lead form slope, etc., could also be considered in this method by adding slicing operation on the gear width direction as described in [26].
As explained in the last paragraph, the first step in building a gear matrix is the determination of the gear tooth profile; three different options could be used to do this: analytically, experimentally or numerically. The analytical method is the most used one in the literature; this method uses the mathematical expression of the involute curve and the tooth root to determine the complete tooth profile curve. The experimental method uses reverse engineering tools and techniques to reconstruct the tooth profile curve from a real gear geometry. Finally, numerical methods use dedicated computer-aided design software to extract the gear tooth profile geometry from the 3D model of the corresponding gear.
The third option is adopted in this study because no explicit mathematical expressions of the teeth root shape are available and experimental methods require sophisticated machines and tools to extract gear teeth profiles from real gears. The CAD models of two commercialized gears were downloaded from the constructor website (KHK-Gears) and the corresponding gear tooth shape was extracted using a CAD software (SOLIDWORKS), as shown in Figure 3. The main differences between the extracted gear tooth geometry and that of theoretical one are illustrated in Figure 3. We can see that the gear addendum diameter is modified, the tooth tip is cornered and the gear tooth width is variable along the tooth length direction. These modifications will affect the gears' contact ratio and the tooth resistance to the applied contact force and, consequently, the TVMS value of the mating gear pair.
The elementary gear tooth profile and side curve are extracted and saved in two vectors, T P and L. An adaptive resampling of these vectors is applied when needed to increase the resolution and the calculation precision, but one should note that excessively high resolutions will lead to high calculation time and so will compromise the method's efficiency. Then, iterative rotational transformations are applied to the T P curve with respect to each tooth's angular position and the resulted vector is assigned to the corresponding tooth column in the gear matrix. Other gear parameters such as addendum and root diameters will also be saved to be used in the next steps.

The Line of Action Determination
In this section, we will discuss the algorithm used to define the line of action equation for arbitrarily chosen coordinates of pinion and gear centers, respectively. The origin of the global frame used in this study is fixed at the pinion reference center point (x r o,p , y r o,p ), so the pinion and gear reference center coordinates are: (x r o,p , y r o,p ) = (0, 0) and (x r o,g , y r o,g ) = (r 1 + r 2 , 0).
As shown in Figure 4, only the gear base circles and centers are required to define the corresponding LOA. Algorithm 1 details the steps followed for the determination of the LOA for varying positions of gear centers: Calculate the pinion-gear center distance: Determine the intersection point of LOA and the line of center distance: where: Figure 4 for more details. 4. Calculate the direction of the LOA: where: α t L is the absolute angle of the LOA in the reference frame of study. Not to be confused with the apparent pressure angle α p as shown in Figure 4.

Contact Point Detection
In this section, we will detail the contact point detection algorithm at each iteration. Let P N×Z be the pinion matrix and G N×Z the gear matrix; this process is realized following these steps:

1.
Set the pinion angular position at time t: where R(θ i ) is the rotation transformation matrix defined as:

2.
Update the position of pinion and gear centers: Determine the line of action equation: Use the algorithm described in Section 2.1.2 to define the line of action equation.

4.
Determine the potential contact points: Sweep the line of action to determine its intersection points with pinion and gear teeth, respectively. Two sets of points are defined at this step: the first set contains the intersection points of the LOA with the pinion teeth and the second one contains the intersection points of the LOA with the gear teeth.
To optimize the exploration time, this operation is limited to the active portion of the line of action. We can prove geometrically that this part is limited by the two addendum circles of the pinion and gear, respectively. 5.
Rearrange the intersection points sets: Use the Euclidean distance to rearrange the two sets of intersection points by crosscomparison; note that the two sets are not necessarily of equal size. This step aims to extract the potential contact point pairs. 6.
Rotate the gear matrix until contact between gear teeth occurs: Use a penalization value ε p to decide whether the points extracted in the previous step are contact points or not (compare the Euclidean distance of potential contact points to ε p ). This penalization value will depend on the resolution used when building the gear matrices and the line of action. If no contact points are detected, the gear matrix is rotated gradually until contact occurs (the penalization constraint is satisfied) as follows: repeat where d s = ±1 depending on the nearest intersection point location and N c is the number of contact points. The gear rotation direction d s is chosen based on the position of the intersection points determined previously. The gear rotation step dθ G should be refined properly to avoid contact jumping, but very small steps could lead to very high computation time. A sufficient condition to avoid this jumping phenomenon is: where r g,h is the gear head radius. 7.
Define contact points: Once contact points at the current iteration are determined, their corresponding indexes, coordinates and the current pressure angle are sent to the mesh stiffness calculation algorithm to determine the corresponding stiffness components. 8.
Update the pinion angular position and return to 1: Increment the pinion angular position and repeat the above described process to determine the contact points at the next iteration.
The algorithm proposed here to locate contact points is not limited to center distance and eccentricity errors only, but it could be used for any positions of pinion and gear centers. Consequently, this algorithm could be used in an online determination of TVMS of gears when simulating the overall dynamic model of the gear system using the positions of pinion and gear centers determined by the system dynamic model integration.

Gear Mesh Stiffness Components
In this section, the outputs of the contact detection algorithm are used in order to calculate the corresponding components of the gear teeth stiffness. In this calculation, the extracted gear tooth geometry is used for the integration of the different stiffness components.
In the potential energy method, the gear tooth is considered as a non-uniform cantilever beam subjected to the contact force F as shown in Figure 5 and the beam theory is applied to calculate the different energies stored in the gear tooth structure. Then, the corresponding stiffnesses are derived from these energies. The stiffness components considered in this study are the Hertzian stiffness k h caused by the local deflection at the contact point position. The tooth bending and shearing stiffnesses (k b , k s ) are caused by the vertical component of the contact force F y . The tooth axial compression stiffness k a is caused by the horizontal component of the contact force F x , and the fillet-foundation stiffness k f is caused by the gear body torsional deformation.
Based on the Hertzian theory, the potential energy method and the work carried out in [27], the mathematical expressions of these components are given as follows: where E, G, ν are the material Young modulus, shear modulus and Poisson coefficient, respectively, and S x , I x are the tooth section area and moment of inertia, respectively, at coordinate x. L f o is the gear foundation width and the coefficients u f , S f , L * , M * , P * , Q * are detailed in [27]. These stiffness elements are evaluated using the discrete trapezoidal integration method. Then, the total mesh stiffness k m at each iteration is determined as follows: where N c is the total number of contact points, i refers to the contact point index and j to the pinion and gear.

Electromechanical Model of a Motor-Gearbox System
To study the effect of assembly errors on the response of a motor-gear system, the electromechanical model of a three-phase induction motor and a one-stage spur gearbox is developed in this section.
A schematic representation of the induction motor is given in Figure 6. This motor is composed of three phases of windings (s a , s b , s c ) on the stator spaced by 120 • and three others on the rotor (r a , r b , r c ). The stator phases are supplied by three-phase sinusoidal voltages at constant frequency and amplitudes, and the rotor windings are short-circuited. The adopted model is based on the following assumptions: proportionality of flux to currents, the influence of skin effect and heating is not taken into account, constant air gap, magneto motor forces are represented with sinusoidal spatial distribution and currents other than in windings are neglected [5]. The application of the fundamental laws of electromagnetic induction gives the following relationships for all windings [28,29]: With R s and R r are the resistance of the stator and rotor winding, respectively. (v sa , v sb , v sc ), (v ra , v rb , v rc ) instantaneous voltages across the stator phases. (i sa , i sb , i sc ), (i ra , i rb , i rc ) instantaneous currents flowing in its phases. (φ sa , φ sb , φ sc ), (φ ra , φ rb , φ rc ) the stator and rotor fluxes.
We observe that each flux interacts with the currents of all the phases including its own. If we take the stator flux on phase "a", the latter can be expressed as follows: ϕ sa = l s i sa + M s i sb + M s i sc + M aa i ra + M ab i rb + M ac i rc (8) where: M sr is the maximum mutual inductance between a phase of the stator and a phase of the rotor. θ m is the angle between the stator frame and the rotor frame. The mutual inductances of the model are dependent on the rotation angle. The Park transformation will be used in order to project the three phases of the windings (a, b, c) of the machine on a frame with two orthogonal two-phase winding (d, q, 0). The purpose of this transformation is to make the mutual inductances independent of the rotation angle. Only the voltage equations on the direct and quadrature axes are used to define the electrical and dynamic model of the induction motor. The equation system constituting the electrical and dynamic model of the induction motor in an equivalent two-phase frame is written as follows: with: ϕ sd = L s i sd + L m i rd ϕ sq = L s i sq + L m i rq and ϕ rd = L r i rd + L m i sd ϕ rq = L r i rq + L m i sq L s is the stator self-inductance. L r is rotor self-inductance. L m is the mutual inductance between stator and rotor. θ s : the electrical angle between the d axis and the stator. θ sl : the electrical angle between the q axis and the rotor. It is therefore possible to estimate from the two preceding equations the stator currents (i sd , i sq ) as well as the rotoric fluxes (φ rd , φ rq ), which will subsequently be used for the estimation of the electromagnetic torque: with ω s = dθ s dt being the stator pulsation and ω sl = ω s − dθ m dt = ω s −θ m , σ = 1 − L m 2 L r L s being the total dispersion coefficient; T r = L r R r : rotor time constant and T s = L s R s : the stator time constant.
The electromagnetic torque can then be expressed as follows: With p being the number of pole pairs. Then, Newton's second law is employed to derive the spur gear dynamic model of Figure 7. For simplification, the LOA direction is assumed to be unchanged in this part of the model, which is justified because we have e i r bi in this study. The reference frame of study is chosen such as the y axis is aligned with the line of action direction as shown in Figure 7, where O i is the rotation center and G i is the geometric center of the gear i.
The equations of motion of the proposed one-stage gearbox are given as follows [9]: where F Ci = m i e iθ 2 i is the centrifugal force and F Ii = m i e iθi is the inertial force as represented in Figure 7. The elastic and damping mesh forces F k and F b are expressed as follows: One can see that the coupling between the mechanical and the electrical models is realized through the motor rotation speedθ m and the electromagnetic torque T em . The proposed model in this study is more precise for simulating the effect of assembly errors because models proposed in the literature often consider only the centrifugal and the inertial efforts induced by the eccentricity and neglect the TVMS deformation or employ approximate formula of this parameter. Additionally, a synchronization is necessary between the centrifugal and inertial fault directions and the TVMS value at each rotation angle; to solve this problem, the TVMS value is injected following the angular positions of the pinion and gear, respectively. This model is general, and one can use the healthy model simply by setting the eccentricities e 1 , e 2 value to zero.  Table 1 shows the geometrical and mechanical parameters of the gear pair used for the TVMS simulation. First, the developed model is used to simulate the gear tooth width modification effect on the final TVMS. Four different configurations were proposed as follows (see Figure 8):
Gear tooth with KHK modification (extracted from the KHK gear 3D model);

3.
Gear tooth with the first proposed modification, a quadratic reduction is proposed with parameters: x m = 2 mm and z m = 1 mm, where x m is the modification length and z m is the modification depth; 4.
Gear tooth with the second proposed modification. a quadratic reduction is proposed with parameters x m = 2 mm and z m = 1.5 mm.
The width modification profiles represented in Figure 8 are the total width variation along the tooth length direction considering the tooth modification from both sides. The TVMS results for the four modifications of the gear tooth width are illustrated in Figure 9. We can see that, globally, the gear tooth width reduction causes a proportional reduction in the gear pair stiffness. For the KHK gear modification, a tiny reduction is observed only in the double-contact period, because the modification length x m covers only the double-contact area of the tooth profile, whereas the two proposed modifications lead to a reduction in the stiffness in the two periods.

Center Distance Error Effect
Center distance is an inevitable error when assembling gear systems. The effect of center distance variation on the final gear mesh stiffness of the studied spur gear pair is investigated in this section. Different stages of the center distance error are considered in this study from d x = 0 mm to d x = 0.8 mm with a step of 0.2 mm. Figure 10 shows the TVMS calculation results for different values of the center distance d x . The simulation results show that increasing the center distance error affects the TVMS quantitatively by reducing the magnitude of both single-contact and double-contact periods and qualitatively by reducing the gears' contact ratio. The variation of the apparent pressure angle and the contact ratio with respect to the center distance variation is detailed in Table 2. One can see that increasing the center distance error causes an increase in the apparent pressure angle and a decrease in the contact ratio. The order spectrum of the TVMS for the different stages of the center distance error is represented in Figure 11. The TVMS signals were calculated in one round. As shown in the TVMS order spectrum, the amplitude of the fundamental mesh frequency decreases when the center distance error increases; however, its harmonics behaves in different ways for each center distance. Indeed, the even harmonics are attenuated for symmetric TVMS curve (the duration of the double-contact phase is equal to that of single-contact phase), and they tend to equal amplitudes when the duration of the double-contact phase becomes very low compared to that of the single-contact phase (Dirac impulse). This means that the center distance error has a direct impact on the system's dynamic response behavior as it is closely related to the TVMS allure.  Figure 12 shows the effect of the center distance error d x in the working part of the gear tooth flank and the distribution of the single-contact and double-contact areas. When the center distance error increases, we can see that the total active area of the tooth flank is reduced, the single-contact area increases and the double-contact area decreases. In addition, one can see that the critical part of the gear tooth moves to the tooth tip direction for greater values of the center distance error; this is very important to note, because almost all methods used for gear strength calculation are based on the location of this critical zone. From the above mentioned results, it was found that the center distance error causes a considerable reduction in the mating gear strength, which will result in higher stresses on the gear teeth. In addition, reducing the contact ratio means that gear teeth will be subjected to these high stresses for extended time periods. Consequently, this will accelerate the degradation and the aging of gear teeth.

Eccentricity Error Effect
The eccentricity error is another common fault occurring when assembling gear systems. It is characterized by a shift between the geometrical center of the gear and that of rotation. This error is incorporated in the previous algorithm by updating the gear center position block using the current gear angular position θ t as follows: x t o,p = −e 1 + e 1 cos(θ t p ) y t o,p = e 1 sin(θ t p ) x t o,g = r 1 + r 2 + e 2 (1 + cos(θ t g )) y t o,g = e 2 sin(θ t g ) These expressions were used to ensure that the center distance would be always greater than r 1 + r 2 , which is the lowest possible value of the center distance. Otherwise, gears will interfere with each other and may not rotate at all. Different configurations of this error are proposed: eccentricities in the pinion, in the gear and in both the pinion and gear and the corresponding impact in the final mesh stiffness are determined. Figure 13 represents the center distance variation with respect to the pinion rotation angle for the three different scenarios considered in this study, an eccentricity of e 1 = 0.2 mm for the pinion and e 2 = 0.25 mm for the gear were considered. As shown in Figures 14 and 15, an eccentricity error causes a double modulation of the TVMS: an amplitude modulation and a frequency modulation. The modulation frequency is equal to that of the defected gear. The amplitude modulation effect is due to the variation of the couple of points being in contact from a mesh period to another due to the modification of the absolute center distance between gears, whereas the frequency modulation is mainly due to the contact ratio modification.   To better illustrate the effect of the eccentricity error in the TVMS, we represented its order spectrum with respect to the pinion rotation angle. To consider the double-eccentricity in the system, the TVMS of the system was calculated in 36 rounds of the pinion which corresponds to a complete period and, for comparison, the remaining configurations were calculated for the same period. Figure 16 shows a close-up of the order spectrum around the fundamental mesh frequency f m = 29. From Figure 16, we can see that an eccentricity error causes a global reduction in the TVMS magnitude proportional to the error severity. In addition, it causes the raise of gears rotation frequencies and sidebands in the order spectrum around the mesh frequency which are equal to the defected gear rotation frequency. These sidebands are caused by the TVMS modulation in amplitude and in frequency.
The computation time of the proposed method for the TVMS calculation depends closely on the discretization rate and the simulation angular step. Two different angular steps are to be fixed: the global resolution step, used to increment the pinion angle from an iteration to another, and a local step, used at each iteration to rotate the gear until contact with pinion occurs. To improve the computation efficiency, the gear is first rotated with a big step proportional to the pinion rotation step (with a coefficient equal to the ratio of the two gear teeth numbers), and then the exploration is executed for contact detection.
Additionally, the gear tooth profile sampling should be refined enough to improve the precision of the discrete integration of the TVMS components. However, the angular step could be adjusted to enhance the computation efficiency; we should note here that this operation is directly related to the dynamic model resolution strategy. If the TVMS is calculated independently and used as an input of the system dynamic model, the calculation resolution could be regulated freely. However, if the TVMS calculation is realized in parallel with the equations of motion integration, the computation time could be very high. To overcome this shortcoming, one can use two different steps during the system resolution-one time step for the equations of motion resolution and another angular step for the TVMS calculation, where the execution of the TVMS model would be conditioned with the pinion angular position regardless of the time step.

Electromechanical Model Simulation Results
The induction motor used for model simulation is a 2.2 kW, 50 Hz, 400 V star-connected, two-pole, 2905 rpm three-phase motor fed directly from the electrical grid. A load torque of 3.7 Nm was applied which corresponds to 50% of the motor full load. The motor and the gearbox parameters used for the system model simulation are detailed in Table 3: The model equations are implemented at Simulink to simulate the electromechanical model response. Figure 17 represents the supply voltage signals, the motor current speed, the developed phase currents and the first pinion displacement y 1 .  Figure 18 represents the first harmonics of the mesh frequency of the displacement signal y 1 spectrum for different stages of the center distance error. We can see that the center distance error value has a direct impact on the total behavior of the displacement spectrum. Consequently, analyzing the form of the signal spectrum could help to detect the presence of the center distance error in the gear system. Figure 19 illustrates the first pinion displacement y 1 for the different configurations of eccentricity error considered in this study. One can see that the presence of an eccentricity error causes a modulation of the displacement signal with the eccentric gear rotation frequency.
Analyzing the displacement signal spectrum represented in Figure 20 allows us to extract the modulation frequencies at each case. As shown in this figure, the eccentricity error leads to a raising of side band peaks around the mesh frequency f m = 1431 Hz and its harmonics, distant with the defected gear frequency f 1 = 49.33 Hz or f 2 = 39.74 Hz depending on which gear is defected (the mesh frequency is theoretically calculated using the rotation frequency and the number of teeth: f m = z 1 f 1 = z 2 f 2 , where f 1,2 are determined from the model resolutionθ 1 ,θ 2 ). In addition, we can see the apparition of the pinion rotation frequency f 1 = 49.33 Hz when it is eccentric because the signal is taken at the pinion position.
To show the eccentricity error effect on the system electrical response, the first phase current spectrum is represented in Figure 21. We can see that in addition to the supply frequency f s = 50 Hz, the mesh frequency f m = 1431 Hz and its harmonics are also modulated by the eccentric gear rotation frequencies f 1 = 49.33 Hz and f 2 = 39.74 Hz.
These simulation results prove the possibility of eccentricity fault detection and isolation by surveying either the vibrational or the electrical response of motor-gearbox systems. In addition, results show that small values of the eccentricity error induce an important deformation of the system's mechanical and electrical responses. These results explain some aspects of the experimental observations of gear vibration and motor phase-currents and prove that assembly errors constitute one of the main reasons for the apparition of side band frequencies around the gear mesh frequencies in addition to other phenomena such as gear tooth profile errors and shaft deformation, etc.   Figure 19. Effect of the eccentricity error on the displacement signal y 1 .
A drawback of these current simulation results is the difficulty to reproduce electrical harmonics on the system current signal; these harmonics could interfere with the fault harmonics, making it difficult to distinguish them in real world measured phase currents, which is mainly due to the adopted model of the induction machine. However, this model enables us to get an overall idea of the effect of these assembly errors on the system electrical response and could contribute a basis to the user when looking for these faults components on the current spectrum.     Figure 21. Effect of the eccentricity error on the phase current signal spectrum.

Conclusions
In this paper, the effect of assembly errors on the gear mesh stiffness of spur gears and the electromechanical response of a motor-gearbox system was determined using a new geometric-based method. A numerical algorithm based on the Euclidean distance minimization was developed to detect actual contact points during the meshing process. This algorithm was then coupled to the potential energy method to calculate the TVMS components corresponding to the extracted contact points. Two types of mounting errors were considered in this study, namely the center distance error and the eccentricity error. The main findings of this study are summarized as follows: • The gear teeth geometry has a direct impact on the gear mesh stiffness curve and modifying the gear teeth width causes a deformation of the TVMS magnitude; • The center distance error considerably affects the TVMS curve, and increasing the center distance causes a global reduction in the TVMS magnitude; it also affects the contact ratio by reducing the proportion of the double-contact period, and consequently, the gear teeth will be subjected to more important stresses; • The effect of single-eccentricity and double-eccentricity errors on the TVMS was determined; it was found that an eccentricity error causes a double modulation of the TVMS amplitude and frequency by the frequency of the eccentric gear; • The center distance error affects the vibration of the system by modulating the magnitudes of the mesh frequency harmonics and the eccentricity error causes the apparition of characteristic side bands around the mesh frequencies separated by the eccentric gear frequency; • In the electrical response of the system, it was found that the presence of an eccentricity error induces a modulation of the mesh frequency components by the corresponding eccentric gear frequency in addition to the modulation by the main supply frequency.
These results could constitute a preliminary basis for detecting and diagnosing this kind of failure, and neglecting these errors could explain some differences between simulation and experiment results in modeling gear dynamics such as the presence of side bands around mesh frequencies. From this perspective, the model of TVMS determination could be further improved to consider gear teeth deformations and load-dependent contact deflections to refine the transition curve between the single-and the double-contact phases. In addition, more advanced models of the induction model could give better results of the effect of these faults on the motor phase current compared to experimental measurements.

Conflicts of Interest:
The authors declare no conflict of interest.

Abbreviations
The following abbreviations are used in this manuscript: