Numerical Modelling of Crack Growth in a Gear Tooth Root

*Corr. Author’s Address: University of Maribor, Faculty of Natural Science and Mathematics, Koroška 160, 2000 Maribor, Slovenia srecko.glodez@uni-mb.si 579 Numerical Modelling of Crack Growth in a Gear Tooth Root Podrug, S. ‒ Glodež, S. ‒ Jelaska, D. Srđan Podrug1 ‒ Srečko Glodež2,* ‒ Damir Jelaska1 1 University of Split, Faculty of Electrical Engineering, Mechanical Engineering and Naval Architecture, Croatia 2 University of Maribor, Faculty of Natural Science and Mathematics, Slovenia


INTRODUCTION
Two kinds of teeth damage can occur on gears under repeated loading due to fatigue; the pitting of gear teeth flanks and tooth breakage in the tooth root [1].In this paper only the tooth breakage is addressed and the developed computational model is used for the calculation of tooth bending strength, i.e. the service life of gear tooth root.
The standardised procedures according to ISO-standards [1] are usually used for an approximate determination of load capacity of gear tooth root.They are commonly based on the comparison of the maximum tooth-root stress with the permissible bending stress.Their determination depends on a number of different coefficients that allow for proper consideration of real working conditions (additional internal and external dynamic forces, contact area of engaging gears, gear material, surface roughness, etc.).The standardised procedures are exclusively based on the experimental testing of the reference gears and they consider only the final stage of the fatigue process in the gear tooth root, i.e. the occurrence of final failure.
However, the complete process of fatigue failure may be divided into the "crack initiation" and "crack propagation" period [2] and [3].An exact definition of the transition from initiation to propagation period is usually not possible.However, the crack initiation period generally accounts for most of the service life, especially in high-cycle fatigue, see Fig. 1.The complete service life of mechanical elements N can than be determined from the number of stress cycles N i required for the fatigue crack initiation and the number of stress cycles N p required for a crack to propagate from the initial to the critical crack length, when the final failure can be expected to occur: The presented work is mainly restricted on the period of the fatigue crack growth, neglecting or merely with experimental determination of the fatigue crack initiation period.In most of recent investigations [3] to [8] a loading cycle of gear meshing is presumed as pulsating acting at the highest point of the single tooth contact.However, in actual gear operation, the magnitude as well as the position of the force, changes as the gear rotates.This fact can be taken into account performing a quasi static numerical simulation in which the gear tooth engagement is broken down into multiple load steps and analyzed separately.In such a way, a more realistic stress cycle in the gear tooth root is obtained resulting in significantly more exact assessment of the crack propagation life, and consequently in the entire fatigue life.In this paper, a similar procedure to the one described in [8] for bevel gears, is used to analyse fatigue crack growth.However, it is appropriately modified and adopted for spur gears.An approach that accounts for fatigue crack closure effects is developed to propagate crack under nonproportional load.

CRACK INITIATION SIZE
In order to calculate the number of stress cycles required for a crack to propagate from the initial to the critical crack length it is necessary to determine fatigue crack initiation size.Although there have been many approaches to determine crack initiation size, there has so far been no perfect approach.One of the most convenient representation of determining the crack initiation size is the Kitagawa-Takahashi plot of applied stress range required for crack growth, Δσ, against crack length, a, using logarithmic scales, as shown in Fig. 2 [9].As the transition point between crack initiation and crack propagation period the threshold crack length a th is selected, below which linear elastic fracture mechanics (LEFM) is not valid.For engineering applications empirical formula for this transition point is proposed [10]: where Δσ FL is the fatigue limit and DK th is the threshold stress intensity range.The threshold crack length a th thus defines the transition point between short and long cracks, i.e. the transition point between initiation and propagation period in engineering applications.However, a wider range of values has been selected for a th in the literature, usually between 0.05 and 1 mm for steels where high strength steels take the smallest values [10] and [11].
a th log a log Dσ Fatigue limit Dσ FL

Short crack regime
Non propagating cracks

Long crack regime (LEFM)
Propagating cracks The application of the linear elastic fracture mechanics (LEFM) to fatigue is based upon the assumption that the fatigue crack growth rate, da/dN, is a function of stress intensity range DK = K max -K min , where a is a crack length and N is a number of load cycles.In this study the simple Paris equation is used to describe the crack growth rate: This equation indicates that the required number of loading cycles N p for a crack to propagate from the initial length a th to the critical crack length a c can be explicitly determined, if C, m and DK(a) are known.C and m are the material parameters and can be obtained experimentally,

Numerical Modelling of Crack Growth in a Gear Tooth Root
usually by means of a three point bending test according to the standard procedure ASTM E 399-80 [12].For simple cases the dependence between the stress intensity range and the crack length DK(a) can be determined analytically as described in [11] and [12].For a more complicated geometry and loading cases it is necessary to use alternative methods.In this work the finite element method (FEM) in the framework of the program package FRANC2D [13], has been used for simulation of the fatigue crack growth, since the uniformly distributed load on the tooth flank is assumed, which enables the usage of two-dimensional finite element mesh.
A different method can be used to determine the equivalent stress intensity range DK eq under mixed mode loading as appears when load is moving on the gear tooth [14] to [22].In presented work the following equation is used to determine the equivalent stress intensity range: where θ 0 is the crack-propagation angle and DK I and DK II are the stress intensity ranges for mode I and mode II, respectively.
To analyse the fatigue crack growth under mix mode conditions the value DK in Eq. (3) has to be replaced with the value DK eq .The crackpropagation angle θ 0 is in this work determined using maximum tensile stress criterion (MTScriterion) as follows: where K I and K II are the stress intensity factors for mode I and mode II, respectively.The complete computational procedure of the fatigue crack propagation under mixed mode loading conditions considering the crack closure effect is described in [3], [23] and [24].

COMPUTATIONAL MODEL
In the proposed computational model, the uniformly load distribution along gear width is assumed, which enables the usage of two-dimensional finite element model.The model is manufactured with the aid of specifically developed software, which on the basis of geometrical parameters determines the rackgenerated gear tooth geometry.In order to capture the correct boundary conditions, one tooth on each side is included in the model.Boundary conditions of the left and right hand edge portions are kept fixed, and since solid gears are explored also the hub portions are kept fixed (Fig. 3).The distance between the root circle and the hub is taken to be of equal tooth height, so that the influence of the fixed hub on tooth base rotation can be neglected.

Fig. 3. FE-model
Two gear models are being explored (Fig. 4): first in which gear tooth is loaded with normal pulsating force acting at the highest point of the single tooth contact (HPSTC), and second in which the fact that in actual gear operation the magnitude as well as the position of the force changes as the gear rotates through the mesh, is taken into account.

Fig. 4. Loading conditions; a) load acting at the HPSTC, b) moving load
The computational analyses are performed for two different gear rim thickness s R in relation to the high of the tooth h (see Fig. 5): (i) s R = 3.3×h, (ii) s R = 0.3×h.

Load Acting at the HPSTC
In that case a loading cycle of meshing gears is presumed as pulsating force acting at the HPSTC (Fig. 4a).The initial crack of length a th calculated using Eq. ( 2) is placed at the critical plane, which is assumed to be perpendicular to the notch surface in a gear tooth root.

Moving Load Model
For a moving load model, a quasi static numerical simulation method is presented in which the gear tooth engagement is broken down into multiple load steps and analyzed separately.During the contact of the teeth pair the load moves along each tooth flank thus changing its direction and intensity.In order to investigate the influence of the moving load on the gear root stress amplitude, the analysis is divided, for example, in sixteen separated load cases (j = 0 to 15) (Fig. 4b).Four of them take the force act on the tooth ahead (0 to 3) and four of them take the force act on the tooth after (12 to 15) the analyzed tooth; in six cases the entire load acts on the analyzed tooth (5 to 10), and in two cases the load is distributed on the two teeth in contact (4 and 11).Force intensity for different load cases can be calculated using the following Eq.: The load sharing factor X Γ which accounts for the load sharing between the various pairs of teeth in mesh along the path of contact for spur gears and no tip relief has a distribution shown in Fig. 6 [1].Γ y is the parameter on the path of contact and can be calculated as follows [1]: where α y is the pressure angle at the treated point Y and α w is the pressure angle at the pitch cylinder.
Fig. 6.Load sharing factor X Γ By analyzing the stress cycle in the gear tooth root it is determined that stress has a maximal value whenever load is in the HPSTC.It follows that the critical plane of the initial crack is a plane perpendicular to the surface at the notch root.The moving load on the gear tooth is nonproportional since the ratio of K II to K I changes during the load cycle.Consequently, the MTScriterion will predict a unique kink angle for each load increment, but in the crack's trajectory is computed at the end of the load cycle.The procedure is fully described in [3].

PRACTICAL EXAMPLE
The crack propagation was analyzed on the gear wheel of the gear pair with basic data given in Table 1.The gear is made of high-strength alloy steel 14CiNiMo13-4 (0.1% C, 0.27% Si, 0.63% Mn, 1.21% Cr, 0.12% Mo, 0.13% Cu, 0.005% P, 0.005% S) with Young's modulus E = 2.07×10 5 MPa, Poison's ratio n = 0.3, ultimate tensile strength R m = 1277 MPa and yield strength R e = 1104 MPa.The gear material is case carburised and ground.Material parameters for crack propagation are given in Table 2.
In numerical computations it has been assumed that the initial crack corresponds to the threshold crack length a th , see Section 1. Considering the material parameters in Table 2 the threshold crack length is equal to a th ≈ 0.02 mm.Its orientation is assumed to be perpendicular to the tooth root surface (Fig. 7).Since crack increment size needs to be prescribed in advance, crack increment size is taken to be 0.005 mm up to the crack length a = 0.2 mm, and after this 0.1 mm to the critical crack length.

Table 1. Basic data of treated spur gear pair
Fig. 8 shows the dependence between the equivalent stress intensity factor K eq and crack length a for two different gear rim thickness (s R = 3.3×h and s R = 0.3×h) if the force is acting at the highest point of the single tooth contact (HPSTC).The FEM-mesh and the crack path for the same cases are shown in Figs. 9 and 10.The similar results are also presented for moving contact loading (see Figs.     16 shows the crack paths which have been determined numerically for different rim thicknesses and different loading conditions.The numerical determined crack paths are then compared with the experimental results taken from [7].A reasonable agreement between numerical and experimental results for deeper rim thickness is observed.This is not the case for thinner rim thickness where the numerical determined crack path significantly differs from the experimental results especially for larger crack lengths.

CONCLUSIONS
The numerical model used to predict the fatigue crack growth in a gear tooth root is Numerical Modelling of Crack Growth in a Gear Tooth Root presented in this paper.The fact that in an actual gear operation the magnitude as well as the position of the force change as the gear rotates through the mesh, is taken into account.In such a way, a more realistic stress cycle in gear tooth root is obtained.The effect of gear rim thickness on the fatigue crack propagation in a gear tooth root and formation of a crack path is also studied.In the numerical computations the crack closure effect is also taken into account, extending an analytical model for plasticity induced crack closure with the partial crack closure concept.In this way, two other closure mechanisms: roughness and oxide induced crack closure are not considered.
Using the numerical procedure described above the predictions of crack propagation lives and crack paths in regard to the gear tooth root stresses are obtained.They are significantly different in comparison to some simplified models, which have been published previously.

Fig. 1 .
Fig. 1.Schematic representation of the service life N of mechanical elements

Fig. 10 .
Fig. 10.The FEM-mesh and crack path for gear rim thickness s R = 0.3×h for HPSTC-loading

Fig. 14 Fig. 11 .Fig. 13 .
Fig.14shows the number of stress cycles N p required for a crack to propagate from the initial (a th ) to the critical (a c ) crack length for gear rim thickness s R = 3.3×h, where two loading conditions are taken into account: (i) normal pulsating force acting at the highest point of the