- Email: [email protected]
Numerical methods for determining stress intensity factors vs crack depth in gear tooth roots
Int..I. Fatigue Vnl, 19, No. I0, pp, ¢~77 685, 1997 *~ 1998 El,,evier Science I,td. All rights reserved Prinled in (;real Britain !)1.4.2 I 123/97/$17.00+.(1II
EI.SEVIER
PII: S0142-1123(97)00101-1
Numerical methods for determining stress intensity factors vs crack depth in gear tooth roots Stanislav Pehan, Trevor K. Hellen*, Jose Flasker and Srecko Glodez University of Maribor, Faculty of Mechanical Engineering, Maribor, Slovenia (Received 7 April 1997; revised 15 July 1997,"accepted 29 July 1997) The calculation of the strength of gears is usually carried out by standard procedures. The calculation of the service life of a gear with a crack in a tooth root. however, is possible only by numerical methods. The first step in determining the service life in such a case is the evaluation of the stress intensity' factor as a function of crack shape and depth. Two-dimensional analysis is appropriate lot this since it is fast and efficient. Here, the direction of crack propagation from the tip of the inilial crack is determined using a special numerical algorithm, whereby the direction of maximum strain energy release rate G is sought. The procedure is repeated incrementally. In order to study the influence of realistic applied loads at the point of contact, the gear has to be treated three-dimensionally. Here. the propagation of each point along the crack tip profile is also assumed to be in the direction of ~he maximum strain energy release rate. The crack depth is determined in such a way' that the stress intensity factor on the crack tip profile is constant. The result of such numerical calculations gives a diagram of the stress intensity factor as a function of crack depth. With known gear material properties it is then possible to calculate the service life of the gear by numerically integrating the Paris equati,m. This article describes two- and three-dimensional methods for monitoring the crack propagation for a particular gear geometry, including the effects of varying through-thickness load behaviour. ~,) 1998 Elsevier Science Ltd. All rights reserved i Keywords: numerical methods; stress intensity factors; gear tooth)
break down. In most cases, this could happen due to the appearance of a fatigue crack in the tooth root. Knowing the process o f crack propagation will help one understand the influence of individual design and technological procedures on the service life of the whole gearbox. As an example of crack propagation, a gear pair of an automotive gearbox installed in at light truck o f 11 000 kg ( 11 tonnes) and 97 kW ( 130 HP) was chosen. Extensive tests were carried out with this vehicle L, so that the loads of the individual gear pair are well known. The basic data of the obserw.'d gear pair are shown in Table 1. The finite element method (FEM) was chosen for the numerical analysis, although some supporting calculations were made using the boundary element method (BEM). The former method was preferred since much more experience was available, and it could deal readily with residual stresses, substructures and fracture p h e n o m e n a as described as follows. BEM has the advantage o f easier modelling since only boundaries and surfaces need be defined. Both FEM two- and three-dimensional analyses are considered, the latter being necessary where the applied load is distributed non-uniformly or the crack is not uniform along the
INTRODUCTION The development in the e c o n o m y o f operation of a transport vehicle means a constant increase in the vehicle payload, higher speed, lower fuel consumption and much longer service life. These requirements also greatly effect the construction of the gearbox. In the 1960s, a truck power of l l 2 k W ( 1 5 0 H P ) was quite normal, but today > 300 kW (400 HP) is required. The gearbox weight is being reduced so that now a gearbox of the same output is twice as light as one in the 1970s. In addition, the expected service life o f the gearbox in the same time has increased from ca 200 000 to > 1 million km. At first sight, the initiation and propagation of a tooth root crack may not seem to be a very significant problem. However, this is probably the only failure mode occurring in such a component. The designer pays great attention to the selection o f gear material, to the technology o f manufacture, and lubricant in order to avoid excessive wear of tooth flanks during the service life. Nevertheless, one day the gear will
*Author for correspondence at: O a k l i e l d House, Thornbury. Bristol B S I 3 1LD. UK.
L o w e r Morton,
677
Stanislav Pehan et al.
678 Table 1
Basic geometric data of pinion and gear
Data Number of teeth Face width Normal module Helix angle Centre distance Working pressure angle Addendum modification
Designation
Pinion
Gear
z~, s2
II 32.5
39 28
b,, b 2 (ram)
m,, (mm) /3 (o) a (mm) c~,, (°)
4.5 () 115 24 0.526
x,, x2
entire width. The translation of the problem to the plane and calculations in two spatial dimensions are considerably quicker and can be automated to a large degree, as appropriate for symmetric loads. For either case, good numerical procedures are described in order to derive sufficiently accurate results. From the resulting curves of stress intensity factor against crack length, Paris' law is used to calculate the number of cycles. Comparisons are included of the two- and three-dimensional results and their implications on the service life. A more detailed description of this last aspect has been published elsewhere 2. LOADS OF FIRST STEP GEARS IN AUTOMOTIVE GEARBOXES The applied load was determined by measuring the load spectrum on the propeller shaft between the gearbox and rear axle differential. The load spectrum for the treated gear pair was determined by multiplying the measured values of the torque with the appropriate gear ratio and also taking into account the share of operation during the entire service life of the vehicle. The stresses calculated according to the equation O"F = O ' F J ~ A K v K F I ~ K F ~ (ISO 6336), are shown in F i g u r e 1, from which it is evident that the face load factor KFt~ has a decisive influence on the stress in the tooth root. The stresses in the tooth root of the gear are always higher than in the pinion root, so consequently only the gear has been analysed. For two-dimensional treatment, the loads are equal to the stresses o-~., whereas in three-dimensions the nominal stresses o'H}
1250,
T
[! ,¢ Gear stress
1000. 750
~,500 250
.......
. . . . . . . . . . . . . . . . . . . . .
factor for gear [
[
.
.
.
.
.
.
.
.
.
.
.
.
.
[]I i, ~
8
0
1 "10 "100 1000 Number of stress repetitions Figure 1
10000
Loads acting on first step gears in an automotive gearbox
0.0593
(KA, Kv and K~;. are equal to 1) present the input information. The influence of the position of the contact area is taken into account by using a real distribution of load there. TWO-DIMENSIONAL TREATMENT OF CRACK PROPAGATION IN A TOOTH ROOT The stress field at the tip of a crack is fully described by the stress intensity factor. Due to the complex geometry of the gear tooth, the stress intensity factor has to be calculated using numerical techniques. Experience in testing the material of the gear with regard to fracture mechanics properties has shown that the gear thickness of 28 mm permits plane strain assumptions. The validity of linear elastic fracture mechanics has to be continually monitored by checking the size of the plastic zone in comparison to the absolute length of the crack. The behaviour of the structure ceases to be linear when the ratio oda exceeds ca 0.4, where ~o is the dimension of plastic zone in the direction of crack propagation and a is the crack length. If this ratio > 0.4, another criterion should be applied 3, since the stress intensity factor loses its validity. Based on the experience of previous research 4, the most suitable package for calculating the crack propagation path, speed and function of stress intensity factor with regard to crack depth is the BERSAFE system 5, based on the FEM. In particular, the virtual crack extension (VCE) method installed in BERSAFE is especially effective. This software was used in conjunction with the IDEAS package, in which useful facilities for the automatic generation of meshes and the colour presentation of stresses, displacements etc. are installed. The nominal load Fp = 1737 N was positioned at the outer point of single engagement. The initial crack was assumed to commence at the point of maximum principal stress on the surface of the tooth root, perpendicular to the surface. The direction of crack propagation at any instantaneous crack tip was determined assuming that the crack propagates in the direction of maximum strain energy release rate. BERSAFE determines this energy corresponding to up to 30 different directions per run, so that the largest value can be determined. The increment of crack propagation length was prescribed initially to ca 0.2 ram, which was gradually increased to 1 mm at step 12, thus defining each successive crack tip. The number of cycles N is governed by Paris' law, which gives daMN = CAK"
( I)
where AK is the stress intensity factor achieved per cycle, and C and n are material constants.
Numerical methods for determining stress intensity factors vs crack depth in gear tooth roots The VCE method calculates first the potential energy release rate, G. This value is then converted into the stress intensity factor. It is assumed that the whole potential energy is converted to the stress intensity factor corresponding to the crack-opening mode I. Mode It is neglected as K , presents only a few percent of the value K~, as was confirmed using BEM~L Gears are case hardened, so that their surfaces have elevated yield stresses in the surface boundary layers. Because of this, the material constants C and n in Equation (1) wiry, and, when case hardened, residual stresses build up in these layers. These properties have been measured for specimens with and without the case hardening treatment, using standard procedures. For an actual gear, however, the main problem is how to model this layer with both sets of material properties present, because it is very thin and contains residual stresses. The residual stresses were simulated using suitable temperatures at each node to represent an equivalent strain field. This simplification gives virtually the same results as would have been obtained by prescribing the precise residual stress field. Two different FEM approaches were used to compare the effects of the residual stresses in the hardened layer. The first approach used material properties (C and n) from a case hardened specimen containing residual stresses distributed exactly as in the real gears. In this case, the residual stresses should not appear in the input to subsequent calculations. In the second approach, C and n are measured from specimens which do not contain residual stresses in their surface layers. Thus, there is no influence of residual stress on the surface layer, and so the residual stress field must be input to the FEM calculations using the above nodal temperature technique. Research showed that the first method did not give correct results, because it was not possible to take into account the influence of residual stresses on the direction of crack propagation. Their contribution in the specimen tests was in determining the material characteristics of the case hardened layer, which only affects the speed of crack propagation. Hence, it is necessary to specify the residual stresses using the second approach. A two-dimensional calculation of crack propagation, taking into account the influence of residual stresses simulated by thermal loads, has been made. The distribution of residual stresses is seen in Figure 2 to predominate in the hardened layer. The initial crack is positioned on the tensile side of the tooth root, perpendicular to the surface. The position of the maximum principal stress is obtained from the combined stress fields of the internal residual stresses, Figure 2, and external load acting on the tooth surface. Figure 3 shows step 12 of the crack propagation from the tooth root. Similar figures were made for every other step of the propagation. The technique for determining the direction of crack propagation is evident from the graph around the top of the crack in Figure 3. This process is repeated until the stress intensity factor, K, reaches its critical value. Especially significant for calculating the service life is a diagram of stress intensity factor as a function of crack length (Figure 4b). The accuracy of the results can be established by comparing the results of different calculation methods. The predicted crack paths, con-
40
80
160
240
679
360
I
Figure 2 Residual stresses due to case hardening and simulated by thermal loading
sisting of the co-ordinates of the individual momentary crack tips, are shown in Figure 4a, obtained from two FEM and one BEM analysis, along with bounding experimental results. Lengths of the plastic ligaments, which appear around the crack tip due to material ductility, indicate whether linear elastic fracture mechanics is valid or not. The relative length of the plastic ligament, compared to the total crack length, rapidly falls from 10% at short cracks, to ca 5%, and just before the final crack length it rises again to ca 8% (Figure 4b). Due to the relatively small average length of plastic ligament, it can be assumed that it does not significantly effect the stress intensity factor, so its influence may be neglected. Otherwise, determining the so-called modified stress intensity factors due to plasticity would require the additional complication of iterative methods. The procedure has thus been simplified and notice taken that errors of the order of 5% in the stress intensity factors may arise from this source. The finite element mesh was designed using local mesh refinements around each instantaneous crack tip, comprising eight triangular finite elements of type EP12F in BERSAFE, each containing the square root displacement singularity. It was also necessary to use two layers of quadrilateral elements EP16 near the surface for forming the case hardened layer. The rest of the mesh was designed so as not to be too dense and elements not too distorted. Triangular elements of type EPI2 were used amongst EPI6 quadrilaterals to aid this. The local refinements around the previous crack tips were left untouched since some effort would have been required without significantly improving the results, although the square root displacement variations were replaced by standard ones. Two different FEM calculations were made, with and without prescribed residual stress fields, as described already. Figure 4 shows a comparison between these two calculations. An additional crack path analysis was performed using a BEM code 6, where the residual stresses were taken into account as a modified material property. There is almost no difference in the calculated stress intensity
Stanislav Pehan et al.
680
Tooth
root / Angle of
4"
crack
propagation (degrees)
r~
_35o •
...... 2,00
Length ofvirtual
-
-
0" 0
°
crack extension 1
1
2
-2"
~
2575
2585
0.01 nun 0.02 mm 0.03 mm 0.04 mm
Figure 3 Crack propagation in the tooth
root, two-dimensional treatment, residual stresses taken into accotlnl. 12th step
#•D..,1737
a b c d
-
limits BEM FEM FEM
of experiments - 2D (two dimensional) - 2D, with residual stresses - 2D, without residual stresses
4000,4.
~
;
,3ooo 4
~\ / C,O
/
"~//
2000,4.1i
/
,.*
~looo
~_~
oi
0.0 ///////////////~v/ P |
~
,is_= .~
,10
Yield .s '-. zone
I.-" "
.20
,
,
.~
2.0
4.0
6.0
Crack length
.
.
.
.
.
.
.
,
o
=
~= "~.
8.0
lmml
b)
Figure4 Results of calculations of crack propagation along the tooth rool, two-dimensional trcatmcm. (a) Crack paths, directions and Ib~ corresponding stress intensity factors
factors but the crack path curves away at quite an early stage. The FEM predictions of crack path are seen to lie within the limiting range of experimental paths, and so seem more plausible. SELECTION OF S U I T A B L E M O D E L FOR THREEDIMENSIONAL TREATMENT In order to obtain an accurate representation of the forces which are exerted on the tooth, including the effects of non-linear distributions down the tooth width, it is necessary to study three-dimensional models. The limited capacity of computer power requires carefully determined Finite Element Models with optimum numbers of elements. With too many elements, the calculation times can become excessive. Two different meshes of the uncracked gear were generated, one with three teeth, Figure 5, and another with one tooth, Figure 6. Different types of elements were compared, including linear and quadratic displacement variations. Although the applied force lies in the two-dimensional plane at the point of single engagement, the actual
load distribution differs down the flank width in the third dimension, so that three different load cases were studied fl)r each set of meshes. Load case A considers a uniform load along the entire tooth width, case B is loaded in the middle part of the width, and case C has a one-sided load near one surface. A comparison of the calculated stresses and displacements of these three load cases is shown in Figures 5 and 6. The displacements of the three teeth FEM models differ from the one-tooth FEM models because the stiffnesses of the two sets of models differ. Consequently, the stress fields would also differ, but probably not so much in the tooth root region. The reference stress field was taken to be that from the model with three teeth using quadratic finite elements. This model would have produced the most accurate results, as confirmed by many other calculations. The Von Mises stresses have been used for basic comparisons. The reference stress field (middle tooth cross section, Figure 7a), is compared to the stress field of the onetooth model (Figure 7b). This comparison shows a good overall agreement, so it was decided to use the
Numerical methods for determining stress intensity factors vs crack depth in gear tooth roots
Load: Fx = 43315 N Fy = 21190 N Solver: S U P E R T A B
Load ease A
LoadcaseB
Load case C Figure 5
- Max. deformation of tooth tip (ram) - Max. tensile stresses in tooth root (NIPs) - Max. Mines' stress in tooth root (MPn) - MnL deformation of tooth tip (ram) - Max. tensile stresses in tooth root (MPa) - Max. Mises' stress in tooth root (MPn) - Max. deformation of tooth tip (main) - Max. tensile stresses in tooth root (MPa) - Max. Mises' stress in tooth root (MPa)
Parabolic Finite EL 4536 Finite El.
21267 nodes
Linear Finite El. 8904 Finite El. 10605 nodes
0.082 0.081 0.080 1126 1176 1169 1017 913 1005 0.102 0.101 0.098 1341 1286 1355 1173 1057 1160 0.120 0.120 0.115 1459 1523 1516 1210 1343 1327
R e s u l t s o b t a i n e d by d i f f e r e n t F E M m e s h e s and d i f f e r e n t e l e m e n t s , t h r e e - d i m e n s i o n a l three-teeth model
Load: Fx = 43315 N Fy = 21190 N
~
A
Load case
B
Load case C
Figure 6
Linear Finite El. 4536 Finite El. 5592 nodes
681
- MaL deformation of tooth tip (mm) - Max. tensile stresses in tooth root (MPa) - MaL Mines' stress in tooth root (MPa) - Max. deformation of tooth tip (ram) - M a L tensile stresses in tooth root (MPa) - Max. MJses' stress in tooth root (MPn) - Max. deformation of tooth tip (ram) - Max. tensile stresses in tooth root (MPa) - Max. Mlses' stress in tooth root (MPa)
r:
~
E
Iver: RSAFE
RTAB
~ s t r s°lver: p ERTAB
Parabolic Finite El. 980 Finite El.
Parabolic Finite El. 980 Finite El.
Parabolic Finite El. 1512 Finite El.
4819 nodes
4819 nodes
7281 nodes
0.067 0.067 0.067 1199 1196 1205 1028 1027 1034 0.085 0.085 0.086 1401 1398 1408 1212 1211 1218 0.102 0.102 0.102 1601 1597 1608 1392 1398 1398
R e s u l t s o b t a i n e d by d i f f e r e n t F E M m e s h e s and d i f f e r e n t e l e m e n t s , t h r e e - d i m e n s i o n a l one-tooth model
model with one tooth and quadratic finite elements l~w all subsequent analyses. Good agreement of the stress fields is particularly evident in the vicinity of the possible crack initiation region, an important consideration when investigating the onset of crack propagation. THREE-DIMENSIONAL TREATMENT OF CRACK PROPAGATION Based on these conclusions, the one-tooth model, which represents the whole gear without any cracks in the tooth root, was used. The main features of the analysis which apply to each of the three load cases are: • Residual stresses are taken into account by specifying suitable temperatures at each node. The resulting thermal stresses represent residual stresses in the case-hardened layer. The field of thermal stresses must be theoretically identical to the field of residual stresses.
@
Restraints are prescribed in such a way that the whole surface of the tooth model which joins on to the gear body is restrained in all directions. In the case of symmetrical load, only one half of the tooth is treated, since full geometrical symmetry exists. Accordingly, restraints are placed on the plane of symmetry by fixing the displacements in the direction of gear thickness (the out-of-plane direction). The tooth model mesh has been divided up using substructure techniques as shown in Figure 8. The one-tooth model is the main structure, whilst the elements around the crack profile constitute the substructure. This division allows much easier manipulation during a sequence of runs when the crack tip profile is being extended, since for each new crack tip position, it is only necessary to change the mesh inside the substructure. The nodes on the boundary between the two structures always remain in the same position. In load case A, the load is uniformly distributed
Stanislav Pehan et al.
682
100
200
300
400
500
600
700
S00
900
1000
MPa
Figure 7
Cross-section in the middle of the tooth: (a) reference model Mises's stresses, load case A, 4536 elements, 21 267 nodes, three teeth: and (b) Mises's stresses, load case A, 980 elements, 4819 nodes, one tooth
Substructure Figure8 Structure clement model
and
substructure
of
the
one
tooth
linitc
along the entire tooth width. It is assumed that the initial crack has a constant depth of 0.3 mm along this width. Each new crack tip profile, a l, a2, etc., in Figure 9, is formed by varying the crack depth through the width to achieve constant, or almost constant, values of stress intensity factor along the profile. This is equivalent to stating that the crack tip profile is correctly modelled when the stress intensity factor is
/-J --
1/'1t I i l i l l
(Tooth center) Figure 9
Calculated fracture surface of the tooth, load case A
constant along it. The modelling of each new crack tip profile involves greater variations of stress intensity factor through the width, and hence greater crack depth variations. This effect is additionally complicated by the need to consider the influence of residual stresses, and because of the plane strain to plane stress transition near the free surface. Figure 10 shows the individual crack tip profiles and corresponding stress intensity factors, and presents the most likely crack propagation under load case A. In this figure and the following two, the gear face axis is along the width, or out-ofplane direction, from 0 to 28ram, with the 14ram point being in the centre of the tooth. For load case B, it is evident from Figure l l b that a reasonably constant value of stress intensity factor is also achieved, therefore the crack tip profile, as shown in Figure 1 la, is probably quite realistic. Figure 12 presents three crack tip profiles fl~r load case C. To obtain a complete picture of crack formation, it is necessary to get additional information about the direction of crack propagation. Such results fin- load case A are shown in Figure 9, which is a threedimensional illustration of the form of the crack tip profile and the direction of crack propagation. A similar illustration fl)r load case C is compared to experimental results in Figure 13, and a satisfactory overall agreement is observed. Finally, Figure 14 shows the dependence of stress intensity factor on the crack depth fbr each point along the tooth thickness for all three load cases. The stress intensity factors obtained by the three-dimensional computations give some scatter, compared to the single curve of the two-dimensional results in Figure 4. Although the latter represents more simplified assumptions, the single curve is much easier for use with the Paris Equation (1) for service life calculations. The three-dimensional scatter is due to varying stress intensity factors along each profile, a realistic effect. [f the central portions of the tooth where the stress intensity factors are constant are considered, then a single curve to which the Paris law could be applied is obtained analogous to the two-dimensional case. By comparing Figures 4 and 14, the two-dimensional results compare well with such a curve from case A, so a similar life prediction is implied. Load case B gives lower stress intensity factors and therefl~re a longer life than the two-dimensional case, whilst load case C is less determinate because less crack depth was achieved. However, the results obtained are higher than the two-dimensional case, indicating a shorter life. To confirm the validity of using linear elastic fracture mechanics, it is necessary to examine the size of plastic zone around the crack tip. An approximate analysis can be made by applying elastic stress fields, but it is better to establish the size of the plastic zone by one of the empirical equations. In the case of plane strain, the size of the plastic zone is calculated from o~ = O.036(K/Rp) 2, and in the plane stress case from ~o = 0.159(K/Rpf. The gear material 20MnCr5 has a yield stress of Rp = 750 MPa and a critical value of stress intensity factor of K,. = 2000 N mm ~P-. K~. is achieved for a crack length of ca 5 mm, as evident from Figure 4. The size of plastic zone determined from these equations amounts to no more than 0.26 ram. This is only ca 5% of the entire crack length, justifying the use of linear elastic fracture mechanics. If the relative
Numerical methods for determining stress intensity factors vs crack depth in gear tooth roots
a) I
I
. . . .
683
b)
a7
i
.
2500
°
2000
~'1500
2
__0
o
:
"
o ~
1000
.~-a4.~-......_~"~
500 14 ~" Figure tO
a l - ~ . ~ - . ~ " 112 IlO 18 6 4 2 Position across gear face (mm)
~
0
0 Position across gear face (mm)
Thrce dimensional load case A results. (a) Crack tips prolilcs, al, a2, elc Ib) ('orrcsponding stre~ ~, intcnsit\ l~tclor
a)
b)
6
1500
..,5 t
1250
41
; ~
a8 ~
2~
~
O
l~---C;~-=-~ 0 i~-------w al ! t
14
12
~' Figure It
a6 a
7,0
5
s00
a2 ~ v
-
i
'10 8 6 14 2 Position across gear face (ram)
!
~ "~ ~
250 0
0
8 Position across gear face (mm)
~
Three-dimensional load case B results. (a) Crack tips proliles, a l , a2, etc. tb~ Corresponding stress, inten'~il\ factor
b) 2500
5-
~-
, 2000
43-
p....-----o-.3
1500
2. 0'-
l
a6
a)
v e-
a8
1000 ,
"=
a7
.~
,ooo
250
1!
0
I
28 -~
24
20 16 12 8 Position across gear face (ram)
4
0 0
[] 28 -~
'24
'20
'16
'12
'8
Position across gear face (ram)
Figure 12 Three-dimensional load case C results. (a) Crack tips profiles, a l , a2, etc. Ibl Corresponding stress intensity factor
size of plastic zone were longer, up to c a 40% of the crack depth, it would also be possible to use linear elastic fracture mechanics. In this case, it would be necessary to make some modifications to account for the extra plasticity. The procedure for determining the size of the yield ligament in three dimensions is similar to the two-dimensional approach, except then this size has to be evaluated along the whole profile. Load case C, where the load is concentrated on one side of tooth width, is considerably more complicated. One-sided loads and one-sided cracks are very fie-
quently found in everyday gear practice. Thus, the reasons permitting the use of a half tooth model in the two previous cases are no longer valid. It is necessary to deal with the entire three-dimensional model of one tooth. Due to software limitations on the semi-bandwidth of the stiffness equations, the number of finite elements in the mesh had to be optimized considerably. The accuracy of such a mesh was confirmed to be adequate by comparing the stress fields of a half tooth model and a whole tooth model. The initial crack was assumed to be as long as the applied
Stanislav Pehan et al.
684
layers, although it was not possible to treat particularly shallow initial cracks using this mesh design.
CONCLUSIONS Analytical methods that are widely standardized enable the calculation of the service life of gears without cracks. Normally, they are sufficient. When a crack in a tooth root appears, the service life of the gear can still be calculated relatively quickly and efficiently for either of the two load cases which are symmetric along the tooth width. The life can be determined using two-dimensional finite element analysis, as described. However, this relatively simple approach fails if it is necessary to consider the influence of a non-symmetric applied load distribution along the tooth root, or the effect of a one-sided crack in a tooth rool with a uniformly distributed load. Load case C is quite common in everyday practice. Then, three-dimensional numerical treatment must be applied as described in this paper. The presented calculations have been shown to be in good agreement with experimental results. Ways of maximizing the efficiency of conducting the relatively complicated three-dimensional analyses ha\e been highlighted. Two particularly useful techniques have been the VCE method, and the use of substructuring for minimizing the amount of mesh through which the crack advances. The same numerical techniques can also be applied to other machine elements its long as a detailed knowledge of the loads and material properties is available.
Figure 13
C a l c u l a t e d tooth fracture s u r f a c e c o m p u r e d with cxpcri ment, load c a s e C
load contact area. As a rule, the crack will propagate in the direction of maximum strain energy release rate. This direction is determined by using the VCE method, as in the two-dimensional case. During each computer run, it is possible to examine 30 different directions of crack extension over any of the nodes lying on the crack profile. The resulting crack depths and stress intensity factors are shown in Figure 12. Since the values of stress intensity factor for each crack tip profile (a l, a2, a3) are reasonably constant along the entire crack width, it could be concluded that the crack tip profiles are meshed correctly. All the crack tip profiles were modelled with equal numbers of element
l
mmn e
2000
•
1500
~
Position across gea r face (case A and B) [ram]:
" ° * ~ A ~"
~
. .'-'A 41~:kO. ,0,
1000
500 0
o
'1
L O A D CASE A
ll0
'3
[] 0.4 00.8
4
'5
6
02.4
2000
A4
1500
A6
×*
Q8
1000 o
500
~
o
,,a
©10 f
o
L O A D CASE B
'1
'3
'4
'5
6
2000 .~
1500 Am
A
1000 500 0
ll0 LOAD CASE C
o
'1
'3
'4
's
Crack depth (mm) F i g u r e 14
×12 .14 Position across gea r face (case C) [mm]:
Slress intensity factors vs c r a c k d e p t h in g e a r toolh roots, load c a s e s A. B a n d C
[] 0.8 O6.1 O 11.4
• 16.6
Numerical methods for determining stress intensity factors vs crack depth in gear tooth roots REFERENCES Tasner, F., Economic Vehicles and Engines, 2nd Part, Research Work. TAM, Maribor, 1985 (in Slovenian). Pehan, S., Hellen, T. K. and Flasker, J., Applying numerical methods for determining the service life of gears. Fatigue Fracture Engineering Material Structures, 1995, 18, 971-979. Schwalbe, K.H., Bruchmechanik Metallischer Werkstg~l'e. Carl Hanser Vcrlag, Munich, 1980.
4 5
6
685
Pehan, S., Velocity of crack propagation in loolh roots considering dynamic loads. MSc. Thesis, University of Maribor, Technical Faculty, Maribor. 1990 (in Slovenian). Hellen, T. K.. On the method of virtual crack extensions. Inter national Joutvud Numerical Methods E,~i,eerin~. 1975, 9, 187-207. Niku, S.M. and Adey, R.A., Crack growth analysis using BE, BENCHmark. NAFEMS .loutvml, March, 1994, 16 19.
EI.SEVIER
PII: S0142-1123(97)00101-1
Numerical methods for determining stress intensity factors vs crack depth in gear tooth roots Stanislav Pehan, Trevor K. Hellen*, Jose Flasker and Srecko Glodez University of Maribor, Faculty of Mechanical Engineering, Maribor, Slovenia (Received 7 April 1997; revised 15 July 1997,"accepted 29 July 1997) The calculation of the strength of gears is usually carried out by standard procedures. The calculation of the service life of a gear with a crack in a tooth root. however, is possible only by numerical methods. The first step in determining the service life in such a case is the evaluation of the stress intensity' factor as a function of crack shape and depth. Two-dimensional analysis is appropriate lot this since it is fast and efficient. Here, the direction of crack propagation from the tip of the inilial crack is determined using a special numerical algorithm, whereby the direction of maximum strain energy release rate G is sought. The procedure is repeated incrementally. In order to study the influence of realistic applied loads at the point of contact, the gear has to be treated three-dimensionally. Here. the propagation of each point along the crack tip profile is also assumed to be in the direction of ~he maximum strain energy release rate. The crack depth is determined in such a way' that the stress intensity factor on the crack tip profile is constant. The result of such numerical calculations gives a diagram of the stress intensity factor as a function of crack depth. With known gear material properties it is then possible to calculate the service life of the gear by numerically integrating the Paris equati,m. This article describes two- and three-dimensional methods for monitoring the crack propagation for a particular gear geometry, including the effects of varying through-thickness load behaviour. ~,) 1998 Elsevier Science Ltd. All rights reserved i Keywords: numerical methods; stress intensity factors; gear tooth)
break down. In most cases, this could happen due to the appearance of a fatigue crack in the tooth root. Knowing the process o f crack propagation will help one understand the influence of individual design and technological procedures on the service life of the whole gearbox. As an example of crack propagation, a gear pair of an automotive gearbox installed in at light truck o f 11 000 kg ( 11 tonnes) and 97 kW ( 130 HP) was chosen. Extensive tests were carried out with this vehicle L, so that the loads of the individual gear pair are well known. The basic data of the obserw.'d gear pair are shown in Table 1. The finite element method (FEM) was chosen for the numerical analysis, although some supporting calculations were made using the boundary element method (BEM). The former method was preferred since much more experience was available, and it could deal readily with residual stresses, substructures and fracture p h e n o m e n a as described as follows. BEM has the advantage o f easier modelling since only boundaries and surfaces need be defined. Both FEM two- and three-dimensional analyses are considered, the latter being necessary where the applied load is distributed non-uniformly or the crack is not uniform along the
INTRODUCTION The development in the e c o n o m y o f operation of a transport vehicle means a constant increase in the vehicle payload, higher speed, lower fuel consumption and much longer service life. These requirements also greatly effect the construction of the gearbox. In the 1960s, a truck power of l l 2 k W ( 1 5 0 H P ) was quite normal, but today > 300 kW (400 HP) is required. The gearbox weight is being reduced so that now a gearbox of the same output is twice as light as one in the 1970s. In addition, the expected service life o f the gearbox in the same time has increased from ca 200 000 to > 1 million km. At first sight, the initiation and propagation of a tooth root crack may not seem to be a very significant problem. However, this is probably the only failure mode occurring in such a component. The designer pays great attention to the selection o f gear material, to the technology o f manufacture, and lubricant in order to avoid excessive wear of tooth flanks during the service life. Nevertheless, one day the gear will
*Author for correspondence at: O a k l i e l d House, Thornbury. Bristol B S I 3 1LD. UK.
L o w e r Morton,
677
Stanislav Pehan et al.
678 Table 1
Basic geometric data of pinion and gear
Data Number of teeth Face width Normal module Helix angle Centre distance Working pressure angle Addendum modification
Designation
Pinion
Gear
z~, s2
II 32.5
39 28
b,, b 2 (ram)
m,, (mm) /3 (o) a (mm) c~,, (°)
4.5 () 115 24 0.526
x,, x2
entire width. The translation of the problem to the plane and calculations in two spatial dimensions are considerably quicker and can be automated to a large degree, as appropriate for symmetric loads. For either case, good numerical procedures are described in order to derive sufficiently accurate results. From the resulting curves of stress intensity factor against crack length, Paris' law is used to calculate the number of cycles. Comparisons are included of the two- and three-dimensional results and their implications on the service life. A more detailed description of this last aspect has been published elsewhere 2. LOADS OF FIRST STEP GEARS IN AUTOMOTIVE GEARBOXES The applied load was determined by measuring the load spectrum on the propeller shaft between the gearbox and rear axle differential. The load spectrum for the treated gear pair was determined by multiplying the measured values of the torque with the appropriate gear ratio and also taking into account the share of operation during the entire service life of the vehicle. The stresses calculated according to the equation O"F = O ' F J ~ A K v K F I ~ K F ~ (ISO 6336), are shown in F i g u r e 1, from which it is evident that the face load factor KFt~ has a decisive influence on the stress in the tooth root. The stresses in the tooth root of the gear are always higher than in the pinion root, so consequently only the gear has been analysed. For two-dimensional treatment, the loads are equal to the stresses o-~., whereas in three-dimensions the nominal stresses o'H}
1250,
T
[! ,¢ Gear stress
1000. 750
~,500 250
.......
. . . . . . . . . . . . . . . . . . . . .
factor for gear [
[
.
.
.
.
.
.
.
.
.
.
.
.
.
[]I i, ~
8
0
1 "10 "100 1000 Number of stress repetitions Figure 1
10000
Loads acting on first step gears in an automotive gearbox
0.0593
(KA, Kv and K~;. are equal to 1) present the input information. The influence of the position of the contact area is taken into account by using a real distribution of load there. TWO-DIMENSIONAL TREATMENT OF CRACK PROPAGATION IN A TOOTH ROOT The stress field at the tip of a crack is fully described by the stress intensity factor. Due to the complex geometry of the gear tooth, the stress intensity factor has to be calculated using numerical techniques. Experience in testing the material of the gear with regard to fracture mechanics properties has shown that the gear thickness of 28 mm permits plane strain assumptions. The validity of linear elastic fracture mechanics has to be continually monitored by checking the size of the plastic zone in comparison to the absolute length of the crack. The behaviour of the structure ceases to be linear when the ratio oda exceeds ca 0.4, where ~o is the dimension of plastic zone in the direction of crack propagation and a is the crack length. If this ratio > 0.4, another criterion should be applied 3, since the stress intensity factor loses its validity. Based on the experience of previous research 4, the most suitable package for calculating the crack propagation path, speed and function of stress intensity factor with regard to crack depth is the BERSAFE system 5, based on the FEM. In particular, the virtual crack extension (VCE) method installed in BERSAFE is especially effective. This software was used in conjunction with the IDEAS package, in which useful facilities for the automatic generation of meshes and the colour presentation of stresses, displacements etc. are installed. The nominal load Fp = 1737 N was positioned at the outer point of single engagement. The initial crack was assumed to commence at the point of maximum principal stress on the surface of the tooth root, perpendicular to the surface. The direction of crack propagation at any instantaneous crack tip was determined assuming that the crack propagates in the direction of maximum strain energy release rate. BERSAFE determines this energy corresponding to up to 30 different directions per run, so that the largest value can be determined. The increment of crack propagation length was prescribed initially to ca 0.2 ram, which was gradually increased to 1 mm at step 12, thus defining each successive crack tip. The number of cycles N is governed by Paris' law, which gives daMN = CAK"
( I)
where AK is the stress intensity factor achieved per cycle, and C and n are material constants.
Numerical methods for determining stress intensity factors vs crack depth in gear tooth roots The VCE method calculates first the potential energy release rate, G. This value is then converted into the stress intensity factor. It is assumed that the whole potential energy is converted to the stress intensity factor corresponding to the crack-opening mode I. Mode It is neglected as K , presents only a few percent of the value K~, as was confirmed using BEM~L Gears are case hardened, so that their surfaces have elevated yield stresses in the surface boundary layers. Because of this, the material constants C and n in Equation (1) wiry, and, when case hardened, residual stresses build up in these layers. These properties have been measured for specimens with and without the case hardening treatment, using standard procedures. For an actual gear, however, the main problem is how to model this layer with both sets of material properties present, because it is very thin and contains residual stresses. The residual stresses were simulated using suitable temperatures at each node to represent an equivalent strain field. This simplification gives virtually the same results as would have been obtained by prescribing the precise residual stress field. Two different FEM approaches were used to compare the effects of the residual stresses in the hardened layer. The first approach used material properties (C and n) from a case hardened specimen containing residual stresses distributed exactly as in the real gears. In this case, the residual stresses should not appear in the input to subsequent calculations. In the second approach, C and n are measured from specimens which do not contain residual stresses in their surface layers. Thus, there is no influence of residual stress on the surface layer, and so the residual stress field must be input to the FEM calculations using the above nodal temperature technique. Research showed that the first method did not give correct results, because it was not possible to take into account the influence of residual stresses on the direction of crack propagation. Their contribution in the specimen tests was in determining the material characteristics of the case hardened layer, which only affects the speed of crack propagation. Hence, it is necessary to specify the residual stresses using the second approach. A two-dimensional calculation of crack propagation, taking into account the influence of residual stresses simulated by thermal loads, has been made. The distribution of residual stresses is seen in Figure 2 to predominate in the hardened layer. The initial crack is positioned on the tensile side of the tooth root, perpendicular to the surface. The position of the maximum principal stress is obtained from the combined stress fields of the internal residual stresses, Figure 2, and external load acting on the tooth surface. Figure 3 shows step 12 of the crack propagation from the tooth root. Similar figures were made for every other step of the propagation. The technique for determining the direction of crack propagation is evident from the graph around the top of the crack in Figure 3. This process is repeated until the stress intensity factor, K, reaches its critical value. Especially significant for calculating the service life is a diagram of stress intensity factor as a function of crack length (Figure 4b). The accuracy of the results can be established by comparing the results of different calculation methods. The predicted crack paths, con-
40
80
160
240
679
360
I
Figure 2 Residual stresses due to case hardening and simulated by thermal loading
sisting of the co-ordinates of the individual momentary crack tips, are shown in Figure 4a, obtained from two FEM and one BEM analysis, along with bounding experimental results. Lengths of the plastic ligaments, which appear around the crack tip due to material ductility, indicate whether linear elastic fracture mechanics is valid or not. The relative length of the plastic ligament, compared to the total crack length, rapidly falls from 10% at short cracks, to ca 5%, and just before the final crack length it rises again to ca 8% (Figure 4b). Due to the relatively small average length of plastic ligament, it can be assumed that it does not significantly effect the stress intensity factor, so its influence may be neglected. Otherwise, determining the so-called modified stress intensity factors due to plasticity would require the additional complication of iterative methods. The procedure has thus been simplified and notice taken that errors of the order of 5% in the stress intensity factors may arise from this source. The finite element mesh was designed using local mesh refinements around each instantaneous crack tip, comprising eight triangular finite elements of type EP12F in BERSAFE, each containing the square root displacement singularity. It was also necessary to use two layers of quadrilateral elements EP16 near the surface for forming the case hardened layer. The rest of the mesh was designed so as not to be too dense and elements not too distorted. Triangular elements of type EPI2 were used amongst EPI6 quadrilaterals to aid this. The local refinements around the previous crack tips were left untouched since some effort would have been required without significantly improving the results, although the square root displacement variations were replaced by standard ones. Two different FEM calculations were made, with and without prescribed residual stress fields, as described already. Figure 4 shows a comparison between these two calculations. An additional crack path analysis was performed using a BEM code 6, where the residual stresses were taken into account as a modified material property. There is almost no difference in the calculated stress intensity
Stanislav Pehan et al.
680
Tooth
root / Angle of
4"
crack
propagation (degrees)
r~
_35o •
...... 2,00
Length ofvirtual
-
-
0" 0
°
crack extension 1
1
2
-2"
~
2575
2585
0.01 nun 0.02 mm 0.03 mm 0.04 mm
Figure 3 Crack propagation in the tooth
root, two-dimensional treatment, residual stresses taken into accotlnl. 12th step
#•D..,1737
a b c d
-
limits BEM FEM FEM
of experiments - 2D (two dimensional) - 2D, with residual stresses - 2D, without residual stresses
4000,4.
~
;
,3ooo 4
~\ / C,O
/
"~//
2000,4.1i
/
,.*
~looo
~_~
oi
0.0 ///////////////~v/ P |
~
,is_= .~
,10
Yield .s '-. zone
I.-" "
.20
,
,
.~
2.0
4.0
6.0
Crack length
.
.
.
.
.
.
.
,
o
=
~= "~.
8.0
lmml
b)
Figure4 Results of calculations of crack propagation along the tooth rool, two-dimensional trcatmcm. (a) Crack paths, directions and Ib~ corresponding stress intensity factors
factors but the crack path curves away at quite an early stage. The FEM predictions of crack path are seen to lie within the limiting range of experimental paths, and so seem more plausible. SELECTION OF S U I T A B L E M O D E L FOR THREEDIMENSIONAL TREATMENT In order to obtain an accurate representation of the forces which are exerted on the tooth, including the effects of non-linear distributions down the tooth width, it is necessary to study three-dimensional models. The limited capacity of computer power requires carefully determined Finite Element Models with optimum numbers of elements. With too many elements, the calculation times can become excessive. Two different meshes of the uncracked gear were generated, one with three teeth, Figure 5, and another with one tooth, Figure 6. Different types of elements were compared, including linear and quadratic displacement variations. Although the applied force lies in the two-dimensional plane at the point of single engagement, the actual
load distribution differs down the flank width in the third dimension, so that three different load cases were studied fl)r each set of meshes. Load case A considers a uniform load along the entire tooth width, case B is loaded in the middle part of the width, and case C has a one-sided load near one surface. A comparison of the calculated stresses and displacements of these three load cases is shown in Figures 5 and 6. The displacements of the three teeth FEM models differ from the one-tooth FEM models because the stiffnesses of the two sets of models differ. Consequently, the stress fields would also differ, but probably not so much in the tooth root region. The reference stress field was taken to be that from the model with three teeth using quadratic finite elements. This model would have produced the most accurate results, as confirmed by many other calculations. The Von Mises stresses have been used for basic comparisons. The reference stress field (middle tooth cross section, Figure 7a), is compared to the stress field of the onetooth model (Figure 7b). This comparison shows a good overall agreement, so it was decided to use the
Numerical methods for determining stress intensity factors vs crack depth in gear tooth roots
Load: Fx = 43315 N Fy = 21190 N Solver: S U P E R T A B
Load ease A
LoadcaseB
Load case C Figure 5
- Max. deformation of tooth tip (ram) - Max. tensile stresses in tooth root (NIPs) - Max. Mines' stress in tooth root (MPn) - MnL deformation of tooth tip (ram) - Max. tensile stresses in tooth root (MPa) - Max. Mises' stress in tooth root (MPn) - Max. deformation of tooth tip (main) - Max. tensile stresses in tooth root (MPa) - Max. Mises' stress in tooth root (MPa)
Parabolic Finite EL 4536 Finite El.
21267 nodes
Linear Finite El. 8904 Finite El. 10605 nodes
0.082 0.081 0.080 1126 1176 1169 1017 913 1005 0.102 0.101 0.098 1341 1286 1355 1173 1057 1160 0.120 0.120 0.115 1459 1523 1516 1210 1343 1327
R e s u l t s o b t a i n e d by d i f f e r e n t F E M m e s h e s and d i f f e r e n t e l e m e n t s , t h r e e - d i m e n s i o n a l three-teeth model
Load: Fx = 43315 N Fy = 21190 N
~
A
Load case
B
Load case C
Figure 6
Linear Finite El. 4536 Finite El. 5592 nodes
681
- MaL deformation of tooth tip (mm) - Max. tensile stresses in tooth root (MPa) - MaL Mines' stress in tooth root (MPa) - Max. deformation of tooth tip (ram) - M a L tensile stresses in tooth root (MPa) - Max. MJses' stress in tooth root (MPn) - Max. deformation of tooth tip (ram) - Max. tensile stresses in tooth root (MPa) - Max. Mlses' stress in tooth root (MPa)
r:
~
E
Iver: RSAFE
RTAB
~ s t r s°lver: p ERTAB
Parabolic Finite El. 980 Finite El.
Parabolic Finite El. 980 Finite El.
Parabolic Finite El. 1512 Finite El.
4819 nodes
4819 nodes
7281 nodes
0.067 0.067 0.067 1199 1196 1205 1028 1027 1034 0.085 0.085 0.086 1401 1398 1408 1212 1211 1218 0.102 0.102 0.102 1601 1597 1608 1392 1398 1398
R e s u l t s o b t a i n e d by d i f f e r e n t F E M m e s h e s and d i f f e r e n t e l e m e n t s , t h r e e - d i m e n s i o n a l one-tooth model
model with one tooth and quadratic finite elements l~w all subsequent analyses. Good agreement of the stress fields is particularly evident in the vicinity of the possible crack initiation region, an important consideration when investigating the onset of crack propagation. THREE-DIMENSIONAL TREATMENT OF CRACK PROPAGATION Based on these conclusions, the one-tooth model, which represents the whole gear without any cracks in the tooth root, was used. The main features of the analysis which apply to each of the three load cases are: • Residual stresses are taken into account by specifying suitable temperatures at each node. The resulting thermal stresses represent residual stresses in the case-hardened layer. The field of thermal stresses must be theoretically identical to the field of residual stresses.
@
Restraints are prescribed in such a way that the whole surface of the tooth model which joins on to the gear body is restrained in all directions. In the case of symmetrical load, only one half of the tooth is treated, since full geometrical symmetry exists. Accordingly, restraints are placed on the plane of symmetry by fixing the displacements in the direction of gear thickness (the out-of-plane direction). The tooth model mesh has been divided up using substructure techniques as shown in Figure 8. The one-tooth model is the main structure, whilst the elements around the crack profile constitute the substructure. This division allows much easier manipulation during a sequence of runs when the crack tip profile is being extended, since for each new crack tip position, it is only necessary to change the mesh inside the substructure. The nodes on the boundary between the two structures always remain in the same position. In load case A, the load is uniformly distributed
Stanislav Pehan et al.
682
100
200
300
400
500
600
700
S00
900
1000
MPa
Figure 7
Cross-section in the middle of the tooth: (a) reference model Mises's stresses, load case A, 4536 elements, 21 267 nodes, three teeth: and (b) Mises's stresses, load case A, 980 elements, 4819 nodes, one tooth
Substructure Figure8 Structure clement model
and
substructure
of
the
one
tooth
linitc
along the entire tooth width. It is assumed that the initial crack has a constant depth of 0.3 mm along this width. Each new crack tip profile, a l, a2, etc., in Figure 9, is formed by varying the crack depth through the width to achieve constant, or almost constant, values of stress intensity factor along the profile. This is equivalent to stating that the crack tip profile is correctly modelled when the stress intensity factor is
/-J --
1/'1t I i l i l l
(Tooth center) Figure 9
Calculated fracture surface of the tooth, load case A
constant along it. The modelling of each new crack tip profile involves greater variations of stress intensity factor through the width, and hence greater crack depth variations. This effect is additionally complicated by the need to consider the influence of residual stresses, and because of the plane strain to plane stress transition near the free surface. Figure 10 shows the individual crack tip profiles and corresponding stress intensity factors, and presents the most likely crack propagation under load case A. In this figure and the following two, the gear face axis is along the width, or out-ofplane direction, from 0 to 28ram, with the 14ram point being in the centre of the tooth. For load case B, it is evident from Figure l l b that a reasonably constant value of stress intensity factor is also achieved, therefore the crack tip profile, as shown in Figure 1 la, is probably quite realistic. Figure 12 presents three crack tip profiles fl~r load case C. To obtain a complete picture of crack formation, it is necessary to get additional information about the direction of crack propagation. Such results fin- load case A are shown in Figure 9, which is a threedimensional illustration of the form of the crack tip profile and the direction of crack propagation. A similar illustration fl)r load case C is compared to experimental results in Figure 13, and a satisfactory overall agreement is observed. Finally, Figure 14 shows the dependence of stress intensity factor on the crack depth fbr each point along the tooth thickness for all three load cases. The stress intensity factors obtained by the three-dimensional computations give some scatter, compared to the single curve of the two-dimensional results in Figure 4. Although the latter represents more simplified assumptions, the single curve is much easier for use with the Paris Equation (1) for service life calculations. The three-dimensional scatter is due to varying stress intensity factors along each profile, a realistic effect. [f the central portions of the tooth where the stress intensity factors are constant are considered, then a single curve to which the Paris law could be applied is obtained analogous to the two-dimensional case. By comparing Figures 4 and 14, the two-dimensional results compare well with such a curve from case A, so a similar life prediction is implied. Load case B gives lower stress intensity factors and therefl~re a longer life than the two-dimensional case, whilst load case C is less determinate because less crack depth was achieved. However, the results obtained are higher than the two-dimensional case, indicating a shorter life. To confirm the validity of using linear elastic fracture mechanics, it is necessary to examine the size of plastic zone around the crack tip. An approximate analysis can be made by applying elastic stress fields, but it is better to establish the size of the plastic zone by one of the empirical equations. In the case of plane strain, the size of the plastic zone is calculated from o~ = O.036(K/Rp) 2, and in the plane stress case from ~o = 0.159(K/Rpf. The gear material 20MnCr5 has a yield stress of Rp = 750 MPa and a critical value of stress intensity factor of K,. = 2000 N mm ~P-. K~. is achieved for a crack length of ca 5 mm, as evident from Figure 4. The size of plastic zone determined from these equations amounts to no more than 0.26 ram. This is only ca 5% of the entire crack length, justifying the use of linear elastic fracture mechanics. If the relative
Numerical methods for determining stress intensity factors vs crack depth in gear tooth roots
a) I
I
. . . .
683
b)
a7
i
.
2500
°
2000
~'1500
2
__0
o
:
"
o ~
1000
.~-a4.~-......_~"~
500 14 ~" Figure tO
a l - ~ . ~ - . ~ " 112 IlO 18 6 4 2 Position across gear face (mm)
~
0
0 Position across gear face (mm)
Thrce dimensional load case A results. (a) Crack tips prolilcs, al, a2, elc Ib) ('orrcsponding stre~ ~, intcnsit\ l~tclor
a)
b)
6
1500
..,5 t
1250
41
; ~
a8 ~
2~
~
O
l~---C;~-=-~ 0 i~-------w al ! t
14
12
~' Figure It
a6 a
7,0
5
s00
a2 ~ v
-
i
'10 8 6 14 2 Position across gear face (ram)
!
~ "~ ~
250 0
0
8 Position across gear face (mm)
~
Three-dimensional load case B results. (a) Crack tips proliles, a l , a2, etc. tb~ Corresponding stress, inten'~il\ factor
b) 2500
5-
~-
, 2000
43-
p....-----o-.3
1500
2. 0'-
l
a6
a)
v e-
a8
1000 ,
"=
a7
.~
,ooo
250
1!
0
I
28 -~
24
20 16 12 8 Position across gear face (ram)
4
0 0
[] 28 -~
'24
'20
'16
'12
'8
Position across gear face (ram)
Figure 12 Three-dimensional load case C results. (a) Crack tips profiles, a l , a2, etc. Ibl Corresponding stress intensity factor
size of plastic zone were longer, up to c a 40% of the crack depth, it would also be possible to use linear elastic fracture mechanics. In this case, it would be necessary to make some modifications to account for the extra plasticity. The procedure for determining the size of the yield ligament in three dimensions is similar to the two-dimensional approach, except then this size has to be evaluated along the whole profile. Load case C, where the load is concentrated on one side of tooth width, is considerably more complicated. One-sided loads and one-sided cracks are very fie-
quently found in everyday gear practice. Thus, the reasons permitting the use of a half tooth model in the two previous cases are no longer valid. It is necessary to deal with the entire three-dimensional model of one tooth. Due to software limitations on the semi-bandwidth of the stiffness equations, the number of finite elements in the mesh had to be optimized considerably. The accuracy of such a mesh was confirmed to be adequate by comparing the stress fields of a half tooth model and a whole tooth model. The initial crack was assumed to be as long as the applied
Stanislav Pehan et al.
684
layers, although it was not possible to treat particularly shallow initial cracks using this mesh design.
CONCLUSIONS Analytical methods that are widely standardized enable the calculation of the service life of gears without cracks. Normally, they are sufficient. When a crack in a tooth root appears, the service life of the gear can still be calculated relatively quickly and efficiently for either of the two load cases which are symmetric along the tooth width. The life can be determined using two-dimensional finite element analysis, as described. However, this relatively simple approach fails if it is necessary to consider the influence of a non-symmetric applied load distribution along the tooth root, or the effect of a one-sided crack in a tooth rool with a uniformly distributed load. Load case C is quite common in everyday practice. Then, three-dimensional numerical treatment must be applied as described in this paper. The presented calculations have been shown to be in good agreement with experimental results. Ways of maximizing the efficiency of conducting the relatively complicated three-dimensional analyses ha\e been highlighted. Two particularly useful techniques have been the VCE method, and the use of substructuring for minimizing the amount of mesh through which the crack advances. The same numerical techniques can also be applied to other machine elements its long as a detailed knowledge of the loads and material properties is available.
Figure 13
C a l c u l a t e d tooth fracture s u r f a c e c o m p u r e d with cxpcri ment, load c a s e C
load contact area. As a rule, the crack will propagate in the direction of maximum strain energy release rate. This direction is determined by using the VCE method, as in the two-dimensional case. During each computer run, it is possible to examine 30 different directions of crack extension over any of the nodes lying on the crack profile. The resulting crack depths and stress intensity factors are shown in Figure 12. Since the values of stress intensity factor for each crack tip profile (a l, a2, a3) are reasonably constant along the entire crack width, it could be concluded that the crack tip profiles are meshed correctly. All the crack tip profiles were modelled with equal numbers of element
l
mmn e
2000
•
1500
~
Position across gea r face (case A and B) [ram]:
" ° * ~ A ~"
~
. .'-'A 41~:kO. ,0,
1000
500 0
o
'1
L O A D CASE A
ll0
'3
[] 0.4 00.8
4
'5
6
02.4
2000
A4
1500
A6
×*
Q8
1000 o
500
~
o
,,a
©10 f
o
L O A D CASE B
'1
'3
'4
'5
6
2000 .~
1500 Am
A
1000 500 0
ll0 LOAD CASE C
o
'1
'3
'4
's
Crack depth (mm) F i g u r e 14
×12 .14 Position across gea r face (case C) [mm]:
Slress intensity factors vs c r a c k d e p t h in g e a r toolh roots, load c a s e s A. B a n d C
[] 0.8 O6.1 O 11.4
• 16.6
Numerical methods for determining stress intensity factors vs crack depth in gear tooth roots REFERENCES Tasner, F., Economic Vehicles and Engines, 2nd Part, Research Work. TAM, Maribor, 1985 (in Slovenian). Pehan, S., Hellen, T. K. and Flasker, J., Applying numerical methods for determining the service life of gears. Fatigue Fracture Engineering Material Structures, 1995, 18, 971-979. Schwalbe, K.H., Bruchmechanik Metallischer Werkstg~l'e. Carl Hanser Vcrlag, Munich, 1980.
4 5
6
685
Pehan, S., Velocity of crack propagation in loolh roots considering dynamic loads. MSc. Thesis, University of Maribor, Technical Faculty, Maribor. 1990 (in Slovenian). Hellen, T. K.. On the method of virtual crack extensions. Inter national Joutvud Numerical Methods E,~i,eerin~. 1975, 9, 187-207. Niku, S.M. and Adey, R.A., Crack growth analysis using BE, BENCHmark. NAFEMS .loutvml, March, 1994, 16 19.