Model-Based Analysis of Spur Gears’ Dynamic Behavior in the Presence of Multiple Cracks

Early detection of tooth cracks is crucial for effective condition-based monitoring and decision making.The scope of this work was to bring more insight into the vibration behavior of spur gears in the presence of single and multiple simultaneous tooth cracks. The investigation was conducted in both time and frequency domains. A finite element analysis was performed to determine the variation in stiffness with respect to the angular position for different combinations of crack lengths. A simplified nonlinear lumped parametermodel of a one-stage gearbox with six degrees of freedomwas then developed to simulate the vibration response of faulty external spur gears. Four differentmultiple-crack scenarios were proposed and studied.The performances of various statistical fault detection indicators were considered and investigated. The simulation results obtained via MATLAB indicated that, as the severity of a single crack increases, the values of the time domain statistical indicators increase also, but at different rates. Moreover, the number of cracks was found to have a negative effect on the values of all the performance indicators, except for the RMS. The number and amplitude of the sidebands in the frequency spectrum were also considered, while assessing the severity of the faults in each scenario. It was observed that, in the case of consecutive tooth cracks, the number of spectrum peaks and the number of cracks were consistent in the frequency range of 4-5 kHz. The main finding of this study was that the peak spectral amplitude was the most sensitive indicator of the number and severity of cracks.


Introduction
Gears are toothed mechanical components that are widely used in numerous industrial applications from heavy machinery to precision instruments to transmit power or motion.In a gear set, regardless of which one is driving the other, the smaller gear is called the pinion, and the larger gear is called the gear or wheel.Gear failure is an alarming and undesirable event that may happen because of an excessive applied load, inadequate lubrication, inaccurate manufacturing, or a bad installation procedure.Gear failure may induce higher unacceptable levels of sound and vibration.It may also decrease the efficiency of transmission, alter the normal operating conditions, and seriously disturb the production rate.In more severe cases, it can also provoke costly consequences that jeopardize machines' safety and even threaten human lives.
Because of more competitive industry conditions, machines are required to work under increasingly extreme operating environments for longer cycles and higher loads.Consequently, the gear teeth become more susceptible to surface fatigue cracks that are usually located in the root of the gear teeth, where the stress is usually at its maximum value.A tooth root crack typically results from insufficient rim thickness in the design, an improperly processed material containing inclusions where cracks can start, or severe operating conditions such as overload or misalignment [1].As shown in Figure 1, several parameters such as the thickness, the width, the length, and the propagation angle () are used to describe a crack.
Condition monitoring of gears is crucially needed to verify their health condition and operational state to detect any potential failure, long time before it becomes a functional failure.Under this maintenance philosophy, usually called condition-based-maintenance, activities are provided by regular periodic data collection, which permits the inspectors to detect degradation before the occurrence of any failure and tells the decision maker, based on the machine's condition, if any maintenance intervention is required.The conditionbased-maintenance approach reduces the risk of failures considerably and improves both the machine's availability and working safety.Typical techniques commonly used in condition monitoring include vibration analysis, oil analysis, particle wear analysis, ultrasonic analysis, infrared thermographic analysis, and motor current signature analysis.Among these methods, vibration-based fault detection and analysis are recognized as the most popular, the most efficient, and the most widely applied in many industries for assessing machine health using the measured vibration signals.Vibration analysis has become highly important in detecting faults in gearing systems.The role of vibration monitoring is to identify any change in the vibration signal caused by gear degradation and to give an early warning.Early gear fault detection allows a proper scheduled shutdown of the whole machine to prevent catastrophic failure [2].However, when traveling from the rotating gear to the sensor, usually installed on the gearbox casing, the vibration signal generated by the meshing of the teeth is always contaminated by vibrations from other different sources.Therefore, collecting the signal requires numerous processing tasks.Vibration signals collected from gearboxes are usually a combination of periodic components (resulting from the interaction of teeth during meshing), transient impacting components (resulting from localized defects on the contacting surfaces), and broadband noise (resulting from the friction between moving surfaces).
External excitations always cause gearbox vibrations from the fluctuation of applied torque and input operating speed and by internal excitations from time-varying cyclic mesh stiffness and transmission errors.The gear mesh stiffness is a significant excitation source for the gear set, so any changes in its value arising from tooth faults can seriously affect the dynamic behavior of the transmission and may lead to an abrupt loss in efficiency.Consequently, numerous investigations have been carried out to analyze and evaluate the mesh stiffness.Most of these studies focus extensively on healthy gears, but for cracked gears, crack modeling and mesh stiffness evaluation are fields that are still currently being explored by many researchers.The finite element method (FEM) and the analytical method (AM) are mainly the two conventional ways utilized for calculating the gear mesh stiffness.The FEM is the most powerful technique used for mesh stiffness evaluation as many gear parameters can be easily incorporated in the model [4].In addition, the results generated by FEM are closer to the real ones.On the other hand, the AM has accomplished a reduction in the computation time and has produced results that agree well with FEM.Analytical investigations of the gear mesh stiffness using crack modeling have been performed by [5] to study the effect of the crack size on the mesh stiffness.The impact of crack propagation in the tooth root on the dynamic response of a gearbox was reviewed in [6].The crack levels were simulated from 0% to 45% of the tooth root thickness.An analytical approach to calculate the reduction in the total gear mesh stiffness caused by the presence of a tooth crack was presented in [3], as well as a model using FEM to verify the results obtained analytically.A modified mathematical model of crack growth in the tooth root was proposed by [7] for calculating time-varying mesh stiffness using an improved potential energy method.In the studies mentioned above, the crack propagation scenario of a constant or uniformly distributed crack length throughout the whole tooth width was assumed.Two additional scenarios have been investigated.The first considered the crack along the whole tooth width with a parabolic length distribution [8]; the second considered crack propagation along both the length and width simultaneously [9].
The dynamic response of gearing systems in the presence of defects (such as spalls or cracks) and gear crack detection can be obtained experimentally using vibration measurement or numerically using simulation models [10].In the literature documenting this research, one can see that a great deal of research has been conducted to analyze the experimentally measured vibration signal for fault detection purposes [11].The advantage of vibration measurements is that they reflect the behavior of a real system, but such measurements are also time consuming and cost consuming, especially when repeated measurements are performed for different crack scenarios.Sometimes, there are limitations in producing real cracks of the right dimensions.
Dynamic modeling and simulation can overcome many of these issues and can be a good alternative for studying the dynamic behavior of a gear system more simply.Dynamic modeling and simulation also have the advantage of increasing the understanding of the system's behavior before the initiation of a measurement campaign.An essential goal of dynamic modeling is to develop a model with a reasonable trade-off between reality and simplicity.During the past few years, gear dynamic modeling has continued to be an alternative approach that is still the subject of much ongoing research.Indeed, different gear dynamic models have been developed and applied for dynamic response simulation, summarized in [12].The lumped mass model is widely used for modeling cracked gear systems such as 4-, 6-, 8-, and 16-DOF models.A gear dynamic model with eight degrees of freedom (DOF) was applied in [13], whereas a six-DOF model was applied in [6], ignoring intertooth friction.A different six-DOF gear dynamic model which considered intertooth friction was conceived in [8].A 16-DOF gear dynamic model was developed by [14] and then adopted by [7] for simulating the dynamic behavior of a one-stage gear system.In addition, FE models also exist and were able to investigate many crack parameters and their effect on the mode shapes, natural frequencies, and the frequency response functions [15].
Although this research area continues to attract the interest of numerous research groups, more focus and research are still needed.Some topics have not yet been covered thoroughly (e.g., issues related to crack modeling, gear mesh stiffness calculation, dynamic modeling, and fault detection methods).The primary objective of the current research is to bring more insights into the understanding of gear dynamics in realistic cases where multiple cracks exist simultaneously in different locations and to varying extents, by using numerical simulation and dynamic modeling.

Research Methodology.
The methodology employed in this study is based on a six-DOF dynamic numerical model [16].It allows the investigation of the effect of one-stage spur gearbox tooth cracks on the vibration response.The contact analysis between the gears was carried out using a tailor-made MATLAB code.The total gear mesh stiffness was estimated with respect to the pinion rotational angle using both SolidWorks and MATLAB software.Tooth root cracks were assumed to be present on the pinion only, with a uniform length extended through the entire tooth width.The total mesh stiffness was then used to simulate the vibration response of the pinion.Intertooth friction is considered in this model.The assumptions used for the development of the dynamic model are similar to those used in [14].It is known that the stiffness of components like bearings and shafts also affects the overall meshing stiffness; however, to reduce the complexity of the proposed model, all the system components, except the gears, were assumed to be rigid, and the influence of lubrication was ignored.The stiffness of the meshing gears was considered, and the error in the mesh stiffness caused by this assumption can be ignored since this study investigates the difference between the healthy and faulty condition.For the vibration analysis, different statistical indicators were applied to the original and residual vibration signals in the time and frequency domains.The diagnostic performance yielded by these statistical techniques (between the original signals and residual signals) was compared and characterized based on their sensitivity.The overall organigram of the code developed is presented in Figure 2.

Contact between Gears
Analysis of the contact between the teeth is an essential step in any dynamic analysis.Knowing how many points are in contact and for how long the teeth are in contact is vital information for determining the gear mesh stiffness.For this calculation, the pinion was supposed to be the driving element.The pinion was considered to be rotating in a counterclockwise direction.As explained in Figure 3, whenever  there is an intersection between a point on the active flank of the tooth and the line of action, this point is entitled to be a contact point.In particular, if that point is located inside the contact region (limits of the meshing zone), it is confirmed as a contact point between the two gears.The graphical results represented in Figure 4, for a gearset with a module of 2 mm and 25 teeth in the pinion and 30 teeth in the gear with a pressure angle of 20 ∘ , clearly show that, during the revolution of both gears, the contact takes place between two teeth at one single point or among four teeth at two different points.It is clear that, during the rotation of the gear, there are two contact points most of the time.The number of contact points and their periodicity depend strictly on the geometry of both mating gears.
Other numbers of teeth and gear ratios were simulated (27 cases in total) with the same module and pressure angle to see their effect on the contact ratio and also to validate the MATLAB code used.One can see that the contact ratio, calculated as the average value of the number of contact points throughout a 360 ∘ rotation, depends not only on the gear ratio (defined as /) but also on the number of teeth on the pinion.The results obtained from the code were compared with the theoretical values calculated using the contact ratio formula [18], and the results were almost identical, as shown in Figure 5, with a percentage error less than 0.2%.

Evaluating Gear Mesh Stiffness with Tooth Cracks
One of the main factors affecting gear mesh stiffness is the crack propagation angle [19].For modeling purposes, the crack is supposed to start from one side at the top of the root fillet and then keeps moving in a straight way towards the center line of the tooth.Once reaching the middle of the tooth, the crack changes its direction and propagates towards the top of the root fillet on the other side.Both propagation lines are straight and have a slope angle of 20 ∘ , as presented in [3,9].Supposing that CL is the length of the crack and PL is the total length of the crack path shown as a red dashed line in Figure 6, the crack length percentage (CLP) can be obtained as The parameters used in the gears modeling are given in Table 1.The backup ratio was taken as 3.3 to avoid the rim thickness effect on the tooth deflection, and root fillet curves were assumed to be circular.
Figure 7 shows the crack propagation path with different CLPs.Within this study, the CLP will be considered to vary from 0% to 45% only.The corresponding crack propagation data are shown in Table 2.
The individual tooth stiffness was evaluated by using a commercial finite element simulation code (SolidWorks), where a "Static" study was performed considering the tooth as a nonuniform beam.Linear-elastic material properties were assumed, as these are reasonable for metal gears.Figure 8  illustrates how the stiffness of a single tooth is calculated, where a force is applied normal to one side of the tooth, acting along the line of action.From the simulation results, the deflections   and   can be used to get  in the direction of the force, which was adopted by [15,16].As the stiffness varies with respect to the angle  between the start of the involute curve and the location at which the force is applied, nine different positions were studied one at a time.The pinion was considered to be fixed on the shaft, as depicted in Figure 9.Each time, the deflection of the tooth is recorded (  ) in the same direction of the force and used to obtain the tooth stiffness at that particular location using the following equation: In Figure 10, the overall mesh of the pinion is displayed.In particular, a finer mesh was used at the location where the force was applied and at the crack region as well.A mesh convergence analysis was conducted, where the mesh element size was decreased until the difference in the deflection value in both  and  directions was almost 2%.The aspect ratio for almost all the elements was less than 3, which avoids any numerical approximation error.The details of the final parameters used are shown in Table 3.
A sample of the final results of the simulation for healthy and faulty pinion is displayed graphically in Figure 11.A healthy gear (no cracks) was simulated as well.The tooth stiffness is plotted against the rotation angle for various crack ratios, both in dimensionless forms: (i) The stiffness ratio [  / max ], where   is the stiffness at position  and  max is the maximum stiffness at the start of the involute curve (bottom of the tooth).(ii) The angle ratio [  / max ], where   is the angle at position  and  max is the total angle of the tooth profile (between the start and the end of the involute curve).
One can see that, when the angle between the contact point and the bottom point (at the start of the involute curve) increases, the distance between the base of the tooth and the point at which the force is applied increases, and consequently the local stiffness at that point decreases.The data points for each case were fitted using a six-degree polynomial curve that approximates the relationship between the stiffness ratio and the angle ratio.The effect of the Hertzian contact ( ℎ ) was taken into consideration as a constant value calculated by the formula presented in [3,16].As previously mentioned, during the meshing between two mating gears, the contact can be single (between two teeth) or double (between two pairs or four teeth) (see Figure 12).Therefore, the total gear mesh stiffness, which is a variable function, could be calculated as follows.For one pair in contact: For two pairs in contact: For the case of a healthy set of gears, the meshing stiffness simulation for a complete cycle of 360 ∘ is displayed in Figure 13.In the graph, one can see the alternation of one pair (low stiffness) and two pairs of teeth being in contact (high stiffness).
To validate the obtained gear mesh stiffness values, a comparison was made with a similar investigation [9] for a healthy and a faulty pinion (one cracked tooth, crack case number (2), length = 0.66 mm, and propagation angle = 20 ∘ ).The results of both studies, portrayed in Figure 14, show a close agreement for both cases.This could be considered as a validation of the approach used in this study.In conclusion, it could be stated that the presence of a gear tooth crack has an adverse effect on the gear mesh stiffness.

Proposed Multiple Cracks Scenarios.
In a practical case, it is very improbable that a single crack would reach 40% Boundary condition "Fixed support around the gear hole"  CLP or more while being the only crack in the entire set of teeth.Usually, when a crack infects one tooth, other cracks are expected to take place on other teeth.These cracks can appear randomly on other teeth, allowing for the possibility of having consecutive and nonconsecutive cracks, as shown in Figure 15.
Different multiple cracks scenarios could be studied, of which the following four scenarios are considered in this work: (i) First scenario: multiple simultaneous cracks with the same length (CLP = 30%) on nonconsecutive teeth (ii) Second scenario: multiple simultaneous cracks with the same length (CLP = 30%) on consecutive teeth (iii) Third scenario: multiple simultaneous cracks with different lengths on nonconsecutive teeth (iv) Fourth scenario: multiple simultaneous cracks with different lengths on consecutive teeth The gear mesh stiffness for each of these previous scenarios is simulated and displayed in Figures 16, 17   respectively.It is clear that there is a difference in the mesh stiffness between consecutive cracks and nonconsecutive cracks.For the case of consecutive cracks, and as expected in advance, the stiffness is lower.

Dynamic Modeling
A six-DOF model was considered in this work, as this is more sensitive to teeth cracks than an 8-or 12-DOF model Total damping between meshing teeth (N s/m) 67 [22] and it was adopted in [8,9,16].The coordinate system was chosen in this model such that one of the axes (the axis) was parallel to the line of action, whereas the -axis was perpendicular to the line of action (see Figure 20).The parameters used in the dynamic model are adopted from [8] and are detailed in Tables 1 and 4. The gears are supported elastically in both directions by springs (1, 2, 1, and 2) and dampers (1, 2, 1, and 2).These elements represent the flexibility introduced by the shafts and the bearings supporting the gears.The radial stiffness and damping of the bearings are considered to be the same both horizontally and vertically.The gearbox casing is supposed to be perfectly rigid.The effect of friction is considered in this study, and the friction coefficient was taken to be 0.06 [8].In fact, the magnitude of the friction force depends on the dynamic friction coefficient () and the contact force (  ) between the teeth along the line of action.The frictional force (  ) is calculated by  However, the transmitted and normal forces have to be determined first.The transmitted force (  ) caused by the torque applied on the pinion (  ) is given as where   is the pitch radius of the pinion.Therefore, the normal force (  ) along the line of action will be obtained via where B is the pressure angle.However, the normal force is not shared equally among the teeth during the varying contact and thus the friction force for the pinion can be calculated as where LSR is the load sharing ratio which varies with the gear rotation and with the crack intensity and was calculated according to the formula presented in [23] as follows: The friction force applied to the gear at the same contact points will be same magnitude as that applied to the pinion but in opposite direction.The friction forces exert moments on the gears.These moments can be calculated, first by identifying the moment arms, taking the geometry of the gear teeth into consideration (Figure 21).
The Cartesian coordinates of the contact  and  are already known from the contact analysis in addition to points  1 2 , so both  1  and  2  can be calculated, and thus the moment arms ( and  for the pinion, and  and  for the gear) can be identified.The angles  1 ,  2 ,  3 , and  4 are calculated using the dot product of two Euclidean vectors as follows:  Thus, the frictional moment arms can be obtained as Finally, the frictional moments are obtained by multiplying the friction force at the contact points with their respective arms.The friction forces and their moments are also considered and explained in [14].
The equations of motion for the system in the -direction are as follows: For the rotary motions of pinion and the gear, the motion equations in the  direction are where  14)-( 19) represent the first and second derivatives with respect to time.For example, ẋ and ẍ represent velocity and acceleration in the  direction, respectively.
It is known that, for some cases of high loading conditions, the tooth meshing stiffness variation may cause variation in operating shaft speed.However, for the sake of simplicity, and because the loading in the present investigation was considered as reasonable, the speed variation of the shafts is neglected.Therefore, the input shaft speed will be precisely fixed at 40 Hz.The numerical solution of the set of equations of motion was achieved with a MATLAB-Simulink code.The main idea of the solving technique was to isolate the term of the higher derivative (acceleration) in the differential equation and to integrate it two times.When the loop is closed, the model will converge to the exact solution after several steps of numerical calculations.The particularity of this method is that it allows reaching a solution even if the system is nonlinear.

Time Domain Analysis.
During the numerical solution of the equations of motion, a fixed-step solver was needed to guarantee obtaining a solution vector with the same length each time.Therefore, different fixed-step ODE Solvers in MATLAB were tested.With a sampling frequency of 400 kHz, all the ODE solvers gave almost the same results except for ODE1, which has a simple integration method [24].The fourth-order Runge-Kutta formula (ODE4) was used for all the simulated cases, as it is widely used for its acceptable accuracy.Normally distributed noise with a signal-to-noise ratio value of 20 dB was added to include the influence of measurement noise [9].In this study, three simulated revolutions were considered.Similarly to any numerical process, a transient response is always encountered before the signal comes to its steady-state regime.Since such transient response is happening mainly in the first revolution, only the second and third periods were kept.The residual signal was then obtained by subtracting the healthy time domain signal from the faulty gear signal to ensure that the remaining signal was only related to the fault.As shown in Table 5, a number of statistical time domain indicators were considered in the present study.Statistical indicators are widely used to evaluate the severity of the fault in a system.By applying these indicators to the obtained signals, faults can be detected at an early stage.Some indicators are traditionally known to have better performance than others, such as the root mean square (RMS) or the kurtosis (KU) [6].In addition, two new parameters, called Talaf and Thikat, previously developed for the condition monitoring of bearings [22], were also considered.
The -displacement of the pinion was analyzed and samples of the time waveform signals, for a healthy case and faulty cases with 25% and 45% CLP, were considered and shown in Figure 22.A noticeable relationship between the crack length and the parameter values can be clearly observed.For the complete topography of the effect of the crack length on the time domain indicators applied to the original signal, the indicators' percentage change was plotted against the CLP (see Figure 23).The indicators' values for a healthy gear were considered as a reference to be used for calculating the indicators' percentage change.Initially, when the crack size was small, the resulting small shocks increased the value of the peak amplitude of the signal, as well as the value of the parameters Thikat, impulse factor (IF), and crest factor (CF).However, they have a minor influence on the RMS, KU, shape factor (SF), and Talaf values.As the size of the crack increases, an increase in the levels of the indicators can be observed.The peak amplitude, Thikat, IF, and CF appear to be more sensitive than the other parameters, whereas the SF appears to be the least sensitive.
Since the values of most of the indicators applied to the original signal did not increase significantly as the CLP increased, the indicators were applied to the residual signal.As expected, the sensitivity of the parameters to the crack length became more pronounced.The percentage change, taking the healthy indicator values as a reference, was calculated and plotted against the CLP, as can be seen in Figure 24.It is apparent that the values of all the statistical parameters increase, as the CLP increased but at different rates.For small crack values, the parameter Thikat had the highest sensitivity.However, when the crack exceeded 25%, kurtosis showed better sensitivity than the other parameters.
These results were compared with another published work [21].However, the crack levels given in [21] are calculated as the crack length divided by the tooth thickness, whereas the crack levels provided in this study are computed as the crack length divided by the entire crack path (3.8 mm).Thus, the crack levels were adjusted to be based on the total crack length.The indicator values were normalized to the healthy signal, and then the kurtosis and crest factor values (applied to the residual signal) obtained by the code were compared with those presented in [21].Figures 25 and 26 indicate a perfect match between the two sets of data point values.
The vibration response of the first and second multiplecrack scenarios was studied to investigate the sensitivity of the parameters above.By plotting the percentage change of the indicators against the number of cracks for the first scenario, the results displayed in Figure 27 clearly demonstrate that, contrary to the previous case, the parameters do not respond the same way when the number of cracks increases.The values of all the parameters, as previously stated, increase because of the presence of one crack.However, as the number of cracks increases, some parameters such as the kurtosis, shape factor, Talaf, and Thikat decrease, and the values for the crest and impact factor will decrease to values even less than that of a healthy signal.On the other hand, the RMS value increases as the number of cracks increases, making it the most sensitive parameter, whereas the peak value remains almost constant after the first crack since the same CLP was used.
In the second scenario, the number of cracked teeth was increased to seven simultaneous cracks, and the values of the statistical parameters were recorded and plotted in Figure 28.All the curves gave the same trends as those of the first scenario.However, the peak amplitude value became constant after three simultaneous cracks.This is because when two cracks happen in consecutive teeth, the vibration response gives higher amplitudes resulting from the dramatic decrease in gear mesh stiffness.
For the third and fourth scenarios, different combinations of crack locations and lengths were considered.From one case to another, the number of cracks and the severity of the existing cracks increased (Case (0) is for a healthy gear).Tables 6 and 7 summarize the details of the simulated cases, where the only difference is the cracks being consecutive or nonconsecutive.The vibration responses of the residual signal for Case (7) for the third and fourth scenarios are shown in Figures 29 and 30, respectively.
If we look at the percentage change of the statistical indicators shown in Figure 31 for the third scenario, it is clear that the peak and Thikat have the highest sensitivity.However, in Cases (8) and above, the value of Thikat starts decreasing.Both CF and the IF are, to some extent, sensitive, but their values decrease after Case (7).A noticeable increase in Talaf and RMS appears after Case (5), where the RMS increases at a higher rate.Kurtosis appears to be more sensitive than RMS but only up to Case (8).
Figure 32 illustrates the indicators' performance for the fourth scenario; almost the same trends are obtained.The main difference is that the parameters Thikat, CF, and IF start decreasing earlier.In general, the sensitivity of all the parameters in this scenario is higher than the previous one.Kurtosis appears to start decreasing after Case (9), which implies that it would start decreasing as well after a few more cracks in the case of the third scenario.The SF, as was concluded before, is the least sensitive parameter and cannot be used to detect the gear tooth cracks.

Frequency Domain Analysis.
The frequency domain analysis is known to have the potential for detecting faults in gears [6,9,23].Other studies have found that the peak   spectral amplitude is more sensitive to gear tooth cracks than the time domain indicators [5,12].Similar to the time domain signal, the signals' spectrum amplitude increases as the severity of the faults increases, and the number of sidebands increases as well [25].Thus, the spectra of all the simulated signals and residual signals were created using two  simulated revolutions for the pinion, for a total of 20,000 samples (10,000 samples per revolution).The spectrum of the healthy case for the original and residual signals is illustrated in Figures 33 and 34, respectively, where they are in close agreement with those presented in [9].Additionally, the gear mesh frequency at 1 kHz and its multiples can be seen in  Since the crack was introduced in the multiples of the rotational speed are expected to be present in the faulty spectrum.Thus, by zooming into Figure 38, harmonics of the pinion rotational speed (40 Hz) can be seen in Figure 39.

RMS
To verify the domain analysis results, the centage change in the maximum peak of the residual signals spectra amplitude for different CLP, obtained from the code, was compared that presented in [9] and it was found to be in good agreement (see Figure 40).
For the proposed multiple-crack scenarios, in the spectra signal of the first scenario (Figure 41), the peaks are sharper than those the scenario (Figure 42) and have more sidebands.Figure 43 shows that as the number of increases, the peak amplitude of the spectra of residual signal for both scenarios linearly, where the rate of the second is higher, as expected.It was also observed that the number of peaks in the residual of the second scenario, between and 5 kHz, corresponds to the number of consecutive cracks in the Figure 42 showed two peaks for two consecutive cracks, and, as illustrated in Figure seven peaks were obtained in the case of consecutive cracks.

Shock and Vibration
Both the number of peaks their amplitudes for the spectra of the residual signal of the third scenario (Figure 45) are higher than that of the fourth scenario (Figure 46).As the number of cracks and their severity increases, the percentage change of the peak of the spectra increases with a third-degree polynomial trend (see Figure 47).

Conclusions
In this study, a numerical model was developed to analyze the dynamic behavior of a one-stage gearbox with external spur gears with an involute tooth profile.A set of MATLAB codes were used to generate the gear tooth profiles, perform the contact analysis, and evaluate the contact ratio.The variable gear mesh stiffness with respect to the angular position was obtained by using the FEM.The total mesh stiffness was then used in a simplified six-DOF nonlinear lumped parameter model to simulate the vibration response of the gears.First, various time domain statistical parameters were extracted from the original, and the residual vibration signals, where the gear was kept healthy, and a single crack was supposed to appear on the pinion.The results of this model were verified in three stages.First, the contact ratios obtained from the contact analysis were compared with the theoretical values, and almost the exact numbers were achieved.Next, the gear mesh stiffness of both a healthy and a cracked pinion was compared with those presented in other research articles, and the results were found to be in good agreement.In the third stage, the values of the statistical indicators applied on both the time and frequency domains at different CLPs were verified using the results of previously published work.
For a more realistic investigation, multiple simultaneous cracks were introduced to the pinion.scenario simulated the effect of the number of multiple cracks up to a total of seven cracks, with the same CLP, located in nonconsecutive teeth.The second scenario is similar to the first scenario, but it considers the possibility of having consecutive cracks.On the other hand, both the third and fourth scenarios study the effect of the number of cracks with different severity, but the fourth scenario simulates consecutive cracks.The gear mesh stiffness of these scenarios  was calculated and inserted in the dynamic model to get the vibration responses.The sensitivity of the previously mentioned indicators was again investigated.
The results show that the parameters have different sensitivity and trends based on the number of cracks, the crack severity, and whether the cracks are consecutive or nonconsecutive.For the first and second scenarios, most of the parameters values were decreasing as the number of cracked teeth increased.However, the RMS value kept increasing as  the number of cracks increased, and the peak value was not significantly affected, as the CLP was constant.If we look at the third and fourth scenarios results, it can be concluded that almost all the indicators increased at first because of the effect of the crack severity but then they start decreasing again as the number of cracks increased further.Contrary to the general trend, the RMS and peak amplitude increased with respect to the growth in severity and number of cracks.It can be observed that the effect of the number of cracks on the statistical indicator parameters is more significant, resulting in a dramatic decrease in the sensitivity of the indicators to crack severity.Therefore, the use of any statistical parameters could be misleading if not considered in the appropriate way.Finally, the peak and the number of sidebands for the frequency domain signal applied to the original and residual signal were investigated.The peak of the residual signal spectra for the first two scenarios was found to increase linearly with the number of cracks.However, the peak value of the third and fourth scenarios was increasing with approximately a third-degree polynomial trend.It was observed that the number of peaks in the residual spectrum of the second scenario between 4 and 5 kHz could be used to predict the number of consecutive cracks.This study has the potential to improve the early detection of gear tooth cracks, as it was found that the spectral amplitude is the most sensitive indicator of the number and severity of cracks.

Figure 1 :
Figure 1: Schematic representation of a crack.(a) A 3D view and (b) a 2D view.

Figure 2 :
Figure 2: The overall organigram for the code developed.

Figure 3 :
Figure 3: Programming flowchart for determining the contact points between two meshing gears.

Figure 5 :
Figure 5: Contact ratio with respect to the gear ratio.

Figure 9 :
Figure 9: Different locations of the force along the tooth surface [17].

Figure 10 :Figure 11 :
Figure10: Overall meshing of the pinion and mesh control at the area where the force was applied and the crack region[17].

Figure 14 :
Figure 14: Comparison of the mesh stiffness for a healthy and a cracked pinion with [9].

Figure 17 :
Figure 17: The gear mesh stiffness for the second scenario (two consecutive cracked teeth with 30% CLP).

Figure 23 :
Figure 23: Performance of different time domain indicators applied to the original signal.

Figure 24 :
Figure 24: Performance of different time domain indicators applied to the residual signal.

Figure 25 :Figure 26 :Figure 27 :Figure 28 :Figure 29 :
Figure 25:  Comparison between the kurtosis values of this study and[21], normalized to the healthy value applied to the residual signal.

Figure 33 .
Figure 33.The spectra of the original signals for a faulty tooth with 25% and 45% crack ratios are shown in Figures 35 and 37, respectively; the residual signals are shown in Figures 36 and 38, respectively.Since the crack was introduced in the multiples of the rotational speed are expected to be present in the faulty spectrum.Thus, by zooming into Figure38, harmonics of the pinion rotational speed (40 Hz) can be seen in Figure39.To verify the domain analysis results, the centage change in the maximum peak of the residual signals spectra amplitude for different CLP, obtained from the code, was compared that presented in[9] and it was found to be in good agreement (see Figure40).For the proposed multiple-crack scenarios, in the spectra signal of the first scenario (Figure41), the peaks are sharper than those the scenario (Figure42) and have more sidebands.Figure43shows that as the number of increases, the peak amplitude of the spectra of residual signal for both scenarios linearly, where the rate of the second is higher, as expected.

Figure 39 :Figure 40 :
Figure 39: Close-up view of Figure 38 showing the multiple integers of the pinion rotational speed (40 Hz).

Figure 41 :Figure 42 :
Figure 41: Spectrum of the residual signal obtained for the first scenario (two nonconsecutive cracked teeth with 30% CLP).

Figure 43 :
Figure 43: Percentage change in the peak of the residual signal spectrum for the first and second scenarios.

Figure 44 :
Figure 44: Spectrum of the residual signal obtained for the second scenario (seven consecutive cracked teeth with 30% CLP).

Figure 45 :Figure 46 :
Figure 45: Spectrum of the residual signal obtained for the third scenario (Case (7) with nonconsecutive cracked teeth).

Figure 47 :
Figure 47: Percentage change in the peak of the residual signal spectra for the third and fourth scenarios.

Table 4 :
[8]ameters of the gear system used in the dynamic model[8].
1 / 2 is mass of the pinion/gear,  1 / 2 is mass moment of inertia of the pinion/gear,  1 / 2 is base circle radius of the pinion/gear,  1 / 2 is horizontal radial stiffness of the input/output bearing,  1 / 2 is vertical radial stiffness of the input/output bearing,  1 / 2 is horizontal radial viscous damping coefficient of the input/output bearing,  1 / 2 is vertical radial viscous damping coefficient of the input/output bearing,  1 / 2 is friction force applied to the pinion/gear,  1 / 2 is friction moment applied to the pinion/gear,  1 is input motor torque,  2 is output torque from load,   is equivalent mesh stiffness,   is mesh damping coefficient,  1 / 2 is linear displacement of the pinion/gear in the -direction,  1 / 2 is linear displacement of the pinion/gear in the -direction, and  1 / 2 is angular displacement of the pinion/gear.Symbols with one or two dots in (

Table 6 :
Cases for the third scenario with different crack locations and lengths.

Table 7 :
Cases for the fourth scenario with different crack locations and lengths.
Four different multiplecrack scenarios were considered in this study.The first