Analytical approach and numerical simulation to investigate the stress field and the dynamic stress intensity factors of a cracked tooth subjected to a periodic loading

A new analytical approach was developed in this paper to study the dynamic analysis of a cracked gear tooth subjected to a periodic loading. The finite elements method (FEM) based on the contour integral technique was used in order to validate the effectiveness and reliability of the analytical formulation. A three-dimensional (3D) model of a cracked tooth was designed where a refined mesh was applied in the crack region to better simulate the stress concentration. The main objective of this study was to investigate the influence of the crack depth and the external load on the variation of the stress intensity factors (SIF) KI and KII associated with the opening crack modes I and II, respectively, and on the stress field near the crack tip. The obtained results show a significant agreement between the analytical results and the FEM findings.


Introduction
Gears, as the most important parts of rotating machinery, are power transmission elements frequently used in industries such as manufacturing, aerospace, automotive, aircraft engines, and so on. They often work in severe conditions and are therefore subject to progressive deterioration especially in the teeth such as wear, cracks, etc. These defects cause a failure in the operation of the mechanical system since they affect the reliability and security of the gear system. Consequently, the study of the dynamic behavior of a cracked gear tooth has been the subject of many investigations owing to its important role in fault detection analysis. Chaari et al. [1] presented an analytical modeling of tooth cracks to quantify the reduction of gear mesh stiffness when the crack is assumed as a straight line shape occurring in the tooth root circle. In order to validate this, a FEM is also computed incorporating the tooth crack. In the same context, Mohammed and Rantatalo [2] developed a dynamic model, including a gyroscopic effect of a one-stage spur gear system, to calculate the time varying mesh stiffness for different crack sizes. Remping et al. [3] investigated the influence of crack position and length on the dynamic characteristics of a cracked gear. Furthermore, the variation of the stress intensity factors related to changes of various parameters was analyzed using the fracture mechanics theory and the finite elements method by ANSYS software. Fakhfakh and Chaari [4] developed a one-stage geared transmission model with two degrees of freedom to study the dynamic response in both cases of healthy and defected teeth. Moreover, they presented a series of tests in the experimental setup to observe and analyze the dynamic behavior and signal responses of perfect and cracked teeth. Nicoletto [5] combined the weight function and the complex potential method in an approximate approach to determine the stress intensity factors of cracked gears taking into account the influence of different parameters such as, loading type, crack propagation direction and teeth number. The modes I and II SIFs were compared with the finite element analysis results reached by reference [6]. An experimental technique based on strain gauge was proposed by Naresh and Anand [7] to measure the gear mesh stiffness of healthy spur gear pair system. Likewise, they calculated the mesh stiffness of cracked spur gear with the help of the strain energy release and the SIF calculation. Thus, the mesh stiffness reduction due to crack have been observed. The experimental results have been validated by the analytical method of Wu et al. [8]. Pehan et al. [9] considered both two-and three-dimensional analysis to evaluate the stress intensity factors in a tooth root as a function of crack depth, and determine the crack propagation direction from the tip of the initial crack.
Zouari et al. [10] presented a finite element method with a three-dimensional survey to study the influence of crack dimension and the direction of crack propagation on the mesh stiffness for spur gear with a crack in the tooth foot. The same computer software (FORTRAN program and ANSYS code) and the same gear geometrical parameters have been used in another work of Zouari et al. [11] in order to estimate the variation of the stress intensity factors according to crack depth and crack propagation angle and track the crack propagation in the tooth foot of a spur gear. Guagliano et al. [12] proposed a versatile approach that considers the contact pressure distribution to evaluate the stress intensity factors along the front of any subsurface cracks found in hypoid gear. Zhanewi et al. [13] established an analytical model based on the potential energy method to investigate the effects on the (TVMS) under three different gear crack types. A three-dimensional (3D) model of spatial cracked gear pair was adopted to verify the reliability of the analytical method. Meanwhile, the effect of these three types of cracks on the equivalent stress, the tooth displacement and contact pressure were analyzed by using the FE method. Recently, Yang et al. [14] proposed analytical formulations, validated by a finite element model, to investigate the influence of the tip fillet radius of rack-type tools and tooth addendum modification on gear mesh stiffness for a spur gear pair.
In light of the review above, the SIF and the stress field around the crack were not extensively studied when a tooth crack appeared and the stress concentration was observed. This needs to be deeply studied. To this end, this paper presented a new analytical approach, solved by applying the Newmark iterative schema in a MATLAB computer program. The proposed approach was validated through a comparison with the FE results provided by a numerical simulation using ABAQUS software. The effects of load and crack depth were discussed, and some conclusions were drawn in the conclusive section.
2 Analytical approach 2.1 Gear tooth modeling A cracked gear tooth, with a width (b) and a crack depth (a), was modeled by a variable section beam fixed at its end. The crack is located at a distance (x c ) from the tooth root as shown in Figure 1.
The load F (t) applied on the tooth is a time varying parameter that reflects gear mesh conditions as the number of teeth in contact varies and as the line of contact of the engaged tooth varies consequently (see Fig. 2). Thus, the tooth is subjected to a periodic load spread over a period (T) and the loading cycle can be modeled by a trapezoidal force (Fig. 3).

Dynamic study of a gear tooth without crack
The tooth was modeled by a recessed beam of a variable section S(x) and a quadratic momentI(x). The dynamic equation of the tooth without crack is written as follows: p(x, t) is a distributed load per unit length, such that: The kinetic energy T and potential strain energy U of the beam can be expressed, respectively, as,   The work achieved by this force F (t) is: The solution to the dynamic equation (2) can be written in the following form: where, Y i ðxÞ¼ x sin ipx 2l is a function which checks the following boundary conditions: Then, we obtain: where, Using the Lagrange equation: we obtain the following equations system: This can be represented in a matrix form as follows, If we limit to i = 1,...., n and j = 1,..., n, we get: qðtÞ ¼ q 1 ðtÞ; :::; q n ðtÞ f g t ; f ðtÞ ¼ f 1 ðtÞ; :::; f n ðtÞ f g t : Applying the time integration schema developed by the Newmark method, and adopting the notation: q ðtÞ ¼ q t , we get: Using the approximation ⃛ q tÀDt ¼ € q t À € q tÀDt Dt we can write: where b ¼ 1 4 and g ¼ 1 2 . Or in the following form: Dt.
The differential equation becomes: It is an equation of the form: where : the fictitious stiffness matrix.

Dynamic analysis of the cracked tooth
In the presence of a crack located orthogonally to the tooth median line, the stress field near the tip of the crack can be expressed in the polar coordinates (r, ') ( Fig. 4) as: K I (t) and K II (t) are the stress intensity factors associated, respectively, with the cracking modes I and II, which can be expressed according to reference [16] by the following relations: where Y I and Y II are factors that depend on the ratio j ¼ a h c between the crack depth and the height of the section. According to [17] we have: where s x (t) and t xy (t) are, respectively, the nominal and the tangential stresses, determined in the case of the healthy tooth and written as, The stress field near the crack tip is given by [18] as shown in the following formulas: t xy ðr; '; tÞ ¼ 3 The finite elements method In order to numerically analyze the dynamic behavior of a cracked tooth and to study the influence of some parameters, namely, the crack depth and the applied load, a three-dimensional model was created and simulated by the ABAQUS software. The cracked tooth was modeled as a non uniform cantilever beam. Assuming that the crack position was set at a distance (x c ) from the tooth root and that this crack extends along the whole tooth width with a uniform depth distribution, the tooth with this type of crack is still considered as a cantilever beam. The boundary conditions of the cracked tooth can be set so that the tooth root is fixed and the top is free. The external load is considered uniformly distributed along the entire tooth width and is located at a distance (x F ) from the tooth root (Fig. 5). The model, for which the geometrical parameters are listed in Table 1, is selected to check the validity of the analytical results.
The finite elements model consists of 17 501 nodes and 15 768 linear hexahedral elements of type C3D8R. To better simulate the stress singularity and obtain accurate results in this paper, the mesh of this model was designed    using local mesh refinements around crack tip and the rest was designed so as not to be too dense and the elements are not too distorted. Figure 6 illustrates the tooth mesh and the refined mesh around crack.
The corresponding steps of the numerical method are illustrated in Figure 7.

Results and discussion
The Contour Integral Method was used in the dynamic implicit simulation to extract the stress field and the stress intensity factors. This method was applied also to clearly show the stress concentration near the crack tip as shown in Figure 8.
The stress intensity factor is a failure criterion that depends on the geometry of the cracked tooth, the size of the crack, as well as the external loading. The dynamic SIF for both mode one and mode two was investigated according to changes in crack length and applied load.   Considering a contact ratio c = 1.6 and rotational speed N = 400 rpm of the pinion which has Z = 15 teeth, the gear mesh period (in second) is defined by [19]: Assuming a crack depth (a = 0.5 mm), while the other parameters are fixed and only the magnitude of the load is varied, the SIFs variation depending on different force values are presented in Figures 9 and 10. Figures 11 and 12 illustrate the SIFs variation for a constant load (F = 500 N), and for various crack depth values with respect to the position of the crack from the tooth root.
From Figures 9-12, it can be observed that the variation of the SIF's (K I and K II ) according to time has a trapezoidal form. This is explained by the periodic force, exerted on the tooth, which varies with time depending on the conditions of the gear and the contact between the teeth.
For the case of low contact ratio (c < 2), a fluctuation of one pair-two pairs of teeth in contact is observed which yields to a time varying stress intensity factor. As it is presented in Figures 9-12, the maximum values of SIF are corresponding to one pair in contact and are observed during (2-c)&T.
The mixed mode of the loading near the crack tip is illustrated with the dynamic SIFs associated with modes I and II. The evolution of stress intensity factors according to the applied load and the crack depth obtained by the numerical method is presented in Figures 13 and 14. It can be observed from these figures that the magnitudes of both    SIFs KI and KII are enhanced with the increasing load and crack depth.
In Figure 14, different crack sizes are selected to study the effect of crack depth on the SIF. As well shown in this figure the SIF's values of both modes I and II were influenced by the growth of crack's size. It also showed that this effect becomes more pronounced when the crack depth exceeds 0.6 mm. Moreover, it can be observed that although KII value increases obviously, it seems always very negligible compared to KI value, and the variation (KI À KII) is still growing.
Furthermore, when a crack appears, it can be seen that the influence of the crack depth on the stress intensity factors is bigger than the effect of the applied load. This is because the cracked tooth loses much rigidity and becomes more and more flexible with the crack depth increase.  Figures 15 and 16 present, under five cases, the effect of crack depth on the stress field near the crack by numerical and analytical methods in order to clearly prove that the stress values will decrease when we are far from the crack. This explains that the trends of the nominal and tangential stresses are parabolic. From these figures, it is clear that the longer the crack is, the more severe the stress concentration becomes. Also, it can be seen that in case 5, which represents a deep crack (a = 1 mm), the levels of s xx and t xy increase significantly compared to that of case 1 which is a smaller crack depth (a = 0.2 mm).
In order to analyze the influence of applied load on the cracked tooth, by two methods, five cases were considered. Each of these cases is related to a different value of force for a constant crack depth (a = 0.5 mm) as shown in Figures 17  and 18. Through Figures 15-18, it can be determined that the magnitudes of the nominal stress (s xx ) and the tangential stress (t xy ) rise with both the load and crack depth increase. Furthermore, the applied load plays a less important role than the crack depth in the growth of stress concentration near the crack. This clearly proves that the stress field changes according to the crack depth rather than to the applied load.
The results of the stress field at the crack tip presented in Figures 15-18 have shown significant agreement between the analytical approach and the numerical findings. These observations demonstrate that this analytical approach will be a major step to study the dynamic behavior for more realistic industrial systems.

Conclusion
Analytical approach and FEM were used to determine the dynamic stress intensity factors for mode I and mode II as well as the stress field in the vicinity of crack, taking into consideration the influence of the crack depth and the variation of applied load. Obviously, the results from the two methods have shown a good agreement and the following conclusions can be forwarded: -The SIF increases with the increase of both the applied load and the crack depth. -The influence of the applied load on the SIF appeared to be smaller than that of the crack depth. -The finite elements analysis, using ABAQUS software and based on the dynamic implicit method, is considered as an efficient and reliable tool to solve dynamic problems and validate the theoretical results. -This analytical approach plays an important role in connecting the fundamental research and the practical usage, and it may be extended in future research studies to describe three-dimensional dynamic problems arising in other types of gearing and defects.