Paper received: 15.09.2009 Paper accepted: 16.05.2011
Numerical Modelling of Crack Growth in a Gear Tooth Root
Srdan Podrug1 - Srecko Glodez2* - 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
A computational model for determination of crack growth in a gear tooth root is presented. Two loading conditions are taken into account: (i) normal pulsating force acting at the highest point of the single tooth contact and (ii) the moving load along the tooth flank. In numerical analysis it is assumed that the crack is initiated at the point of the largest stresses in a gear tooth root. The simple Paris equation is then used for a further simulation of the fatigue crack growth. The functional relationship between the the stress intensity factor and crack length K = f(a), which is needed for determining the required number of loading cycles N for a crack propagation from the initial to the critical length, is obtained using a displacement correlation method in the framework of the FEM-method considering the effect of crack closure. The model is used for determining fatigue crack growth in a real gear made from case carburised and ground steel 14CiNiMo13-4, where the required material parameters were determined previously by appropriate test specimens. The results of the numerical analysis show that the prediction of crack propagation live and crack path in a gear tooth root are significantly different for both loading conditions considered.
© 2011 Journal of Mechanical Engineering. All rights reserved. Keywords: gears, fatigue, crack growth, numerical modelling
0 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 Ni required for the fatigue crack initiation and the number of stress cycles Np required for a crack to propagate from the initial to the critical crack length, when the final failure can be expected to occur:
N = Ni + Np.	(1)
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
*Corr. Author's Address: University of Maribor, Faculty of Natural Science and Mathematics, Koroška 160, 2000 Maribor, Slovenia srecko.glodez@uni-mb.si
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.
Crack propagation period Vp
log N
Fig. 1. Schematic representation of the service life N of mechanical elements
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.
not valid. For engineering applications empirical formula for this transition point is proposed [10]:
1
ath
AK«
Act
FL
(2)
where AcFL is the fatigue limit and DKth is the threshold stress intensity range. The threshold crack length ath 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 ath 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
Fatigue limit AaFL
Propagating cracks
Non propagating cracks
Short crack regime
Long crack regime (LEFM)
ath	log a
Fig. 2. Kitagawa-Takahashi plot
2 FATIGUE CRACK PROPAGATION
1 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, Ac, 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 ath is selected, below which linear elastic fracture mechanics (LEFM) is
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 = Kmax-Kmin, 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:
— = C [AK (a)]" d N L J
(3)
This equation indicates that the required number of loading cycles Np for a crack to propagate from the initial length ath to the critical crack length ac can be explicitly determined, if C, m and DK(a) are known. C and m are the material parameters and can be obtained experimentally,
2
usually by means of a three point bending test according to the standard procedure ASTM E 39980 [12]. For simple cases the dependence between the stress intensity range and the crack length AK(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 AKeq 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:
AKeq = cos2 ——
eq 2
AK T cos—^ - 3AKTTsin-^-
1 2 TT 2
(4)
where 60 is the crack-propagation angle and AKI and AKII are the stress intensity ranges for mode I and mode II, respectively.
To analyse the fatigue crack growth under mix mode conditions the value AK in Eq. (3) has to be replaced with the value AKeq. The crack-propagation angle 80 is in this work determined using maximum tensile stress criterion (MTS-criterion) as follows:
= 2 tan
1. Kl
4 K
f 2
{ K TT y
+ 8
(5)
where KT and KTT 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].
3 COMPUTATIONAL MODEL
Tn 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 rack-generated gear tooth geometry. Tn 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 sR in relation to the high of the tooth h (see Fig. 5):
(i)	sR = 3.3*h,
(ii)	sR = 0.3xh.
Fig. 5. Different rim thickness of analysed gear; a) sR = 3.3*h, b) sR = 0.3*h
3.1	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 ath 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.
3.2	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.:
Fj - FHpSTC ■ XT.
(6)
The load sharing factor Xr 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]. r is the parameter on the path of contact and can be calculated as follows [1]:
tana,,
r
tan a
--1,
(7)
where ay is the pressure angle at the treated point Y and aw is the pressure angle at the pitch cylinder.
Fig. 6. Load sharing factor Xr
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 Kn to Ki changes during the load cycle. Consequently, the MTS-criterion 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].
4 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.07x105 MPa, Poison's ratio v = 0.3, ultimate tensile strength Rm = 1277 MPa and yield strength Re = 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 ath, see Section 1. Considering the material parameters in Table 2 the threshold crack length is equal to ath ~ 0.02 mm.
Its orientation is assumed to be perpendicular to the tooth root surface (Fig. 7).
Table 1. Basic data of treated spur gear pair [3]
Magnitude	Value
Number of teeth for pinion	z1 = 28
Number of teeth for wheel	z2 = 28
Module	mn = 3.175 mm
Addendum modification	xj = -0.05
coefficient for pinion	
Addendum modification	x2 = -0.05
coefficient for wheel	
Gear width for pinion	b1 = 6.35 mm
Gear width for wheel	b2 = 6.35 mm
Flank angle of tool	an = 20°
Radial clearance factor	c* = 0.35
Relative radius of curvature	pf* = 0.35
of tool tooth	
Addendum of tool	ha* = 1.05
Dedendum of tool	h* = 1.35
Tip diameter	Standard clearance
Table 2. Material parameters for crack propagation [3]
Magnitude	Value
Threshold stress intensity range	AKth = 122 Nmm-3/2
Fracture toughness	KIc = 2954 Nmm-3'2
Material parameter of Paris equation	C = 3.128x10-13
Material exponent of Paris equation	m = 2.954
Fatigue limit	Agfl = 450 MPa
Fig. 7. Initial crack orientation
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.
Fig. 8 shows the dependence between the equivalent stress intensity factor Keq and crack length a for two different gear rim thickness (sR = 3.3*h and sR = 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. 11 to 13).
	3500
	3000
°E	2500
nl	
Oh	
2	2000
	
	1500
	1000
	500
											
	ïK - 3 3 vfc ___sR = 0.3.-A										
											
											
											
											
											
											
											
0 0.5 1.0 1.5 2.0 2.5 3.0 3.5 4.0
4.5 5.0 rmml
Fig. 8. The diagram (Keq - a) for HPSTC-loading
Fig. 9. The FEM-mesh and crack path for gear rim thickness sr = 3.3*h for HPSTC-loading
Fig. 10. The FEM-mesh and crack path for gear rim thickness sr = 0.3*h for HPSTC-loading
Fig. 14 shows the number of stress cycles Np required for a crack to propagate from the initial (ath) to the critical (ac) crack length for gear rim thickness sR = 3.3*h, where two loading conditions are taken into account: (i) normal pulsating force acting at the highest point of the
single tooth contact (HPSTC), and (ii) the load moves along the tooth flank. It is clear that the crack grows faster in the case of moving loading conditions. Similar results for gear rim thickness
4.0 4.5 5.0 a [mm]
Fig. 11. The diagram (Keq - a) for moving contact loading
Fig. 12. The FEM-mesh and crack path for gear rim thickness Sr = 3.3xh for moving contact loading
Fig. 13. The FEM-mesh and crack path for gear rim thickness Sr = 0.3xh for moving contact loading
Fig. 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.
10 12 14 16 N-104 [cycles]
Fig. 14. Crack propagation live for Sr = 3.3xh
ÀM04 [cycles] Fig. 15. Crack propagation live for sR = 0.3*h
Fig. 16. Crack paths for; a) Sr = 3.3xh, b) SR=0.3*h (A = numerical for HPSTC loading, B = numerical for moving loading, C = experimental)
5 CONCLUSIONS
The numerical model used to predict the fatigue crack growth in a gear tooth root is
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.
6 REFERENCES
[1]	ISO 6336 (2006). Calculation of load capacity of spur and helical gears, International Standard, Geneve.
[2]	Glodez, S., Sraml, M., Kramberger, J. (2002). A computational model for determination of service life of gears. International Journal of Fatigue, vol. 24, p. 1013-1020.
[3]	Podrug, S., Jelaska D., Glodez, S. (2008). Influence of different load models on gear crack path shapes and fatigue lives. Fatigue and Fracture of Engineering Materials and Structures, vol. 31, p. 327-339.
[4]	Pehan, S., Hellen, T.K., Flasker, J., Glodez, S. (1997). Numerical methods for determining stress intensity factors vs crack depth in gear tooth root. International Journal of Fatigue, vol. 19, p. 677-685.
[5]	Blarasin, A., Guagliano, M., Vergani, L. (1997). Fatigue crack growth prediction in specimens similar to spur gear teeth. Fatigue and Fracture of Engineering Materials and Structures, vol. 20, p. 1171-1182.
[6]	Kato, M., Deng, G., Inoue, K., Takatsu, N. (1993). Evaluation of the strength of carburized spur gear teeth based on fracture mechanics. JSME International Journal, vol. 36, p. 233-240.
[7]	Lewicki, D.G., Ballarini, R. (1997). Rim thickness effects on gear crack propagation life. International Journal of Fatigue, vol. 87, p. 59-86.
[8]	Spievak, L.E., Wawrzynek, P.A., Ingraffea, A.R., Lewicki, D.G. (2001). Simulating fatigue crack growth in spiral bevel gears.
Engineering Fracture Mechanics, vol. 68, p. 53-76.
[9]	Kitagawa, H., Takahashi, S. (1976). Applicability of fracture mechanics to very small cracks or cracks in the early stage. Proceedings of the 2nd International Conference on the Behaviour of Materials, p. 627-631.
[10]	Bhattacharya, B., Ellingwood B. (1998). Continuum damage mechanics analysis of fatigue crack initiation. International Journal of Fatigue, vol. 20, p. 631-639.
[11]	Ewalds, H.L., Wanhill, R.J. (1989). Fracture Mechanics. Edward Arnold Publication, London.
[12]	ASTM E 399-80, American standard, West Conshohocken.
[13]	FRANC2D (2000). User's Guide, Version 2.7. Cornell University, Ithaca
[14]	Shih, C.F., de Lorenzi, H.G., German, M.D. (1976). Crack extension modelling with singular quadratic isoparametric elements. International Journal of Fracture, vol. 12, p. 647-651.
[15]	Narayana, K., Dattaguru, B. (1996). Certain aspects related to computation by modified crack closure integral (MCCI). Engineering Fracture Mechanics, vol. 55, p. 335-339.
[16]	Raju, I.S., Shivakumar, K.N. (1990). An equivalent domain integral method in the two-dimensional analysis of mixed mode crack problems. Engineering Fracture Mechanics, vol. 37, p. 707-725.
[17]	Bittencourt, T.N., Wawrzynek, P. A., Ingraffea, A.R., Sousa, J.L. (1996). Quasi-automatic simulation of crack propagation for 2D LEFM problems. Engineering Fracture Mechanics, vol. 55, p. 321-334.
[18]	Erdogan, F., Sih, G.C. (1963). On the crack extension in plates under plane loading and transverse shear. Journal of Basic Engineering D, vol. 85, p. 519-525.
[19]	Sih, G.C. (1974). Strain energy density factor applied to mixed mode crack problems.
International Journal of Fracture, vol. 10, p. 305-321.
[20]	Hussain, M.A., Pu, S.L., Underwood, J. (1974). Strain energy release rate for a crack under combined mode I and mode II. Fract Anal ASTMSTP. vol. 560, p. 2-28.
[21]	Yan, X., Du, S., Zhang, Z. (1992). Mixed-mode fatigue crack growth prediction in biaxially streched sheets. Engineering Fracture Mechanics, vol. 43, p. 471-475.
[22]	Abdel Mageed, A.M., Pandey, R.K. (1992). Studies on cyclic crack path and the mixed-mode crack closure behaviour in Al alloy.
International Journal of Fatigue, vol. 14, p. 21-29.
[23]	Budiansky, B., Hutchinson, J.W. (1978). Analysis of closure in fatigue crack growth.
Journal of Applied Mechanics, vol. 45, p. 267-276.
[24]	Kujawski, D. (2001). Enhanced model of partial crack closure for correlation of R -ratio effects in aluminum alloys. International Journal of Fatigue, vol. 23, p. 95-102.