Rolling contact fatigue crack path prediction by the asperity point load mechanism

Engineering Fracture Mechanics 78 (2011) 2848–2869

Contents lists available at SciVerse ScienceDirect

Engineering Fracture Mechanics journal homepage: www.elsevier.com/locate/engfracmech

Rolling contact fatigue crack path prediction by the asperity point load mechanism D. Hannes, B. Alfredsson ⇑ Department of Solid Mechanics, KTH Royal Institute of Technology, Stockholm, Sweden

a r t i c l e

i n f o

Article history: Received 17 December 2010 Received in revised form 19 May 2011 Accepted 29 July 2011

Keywords: Rolling contact fatigue Spalling Asperity Fatigue crack path Mode I Plane mixed-mode

a b s t r a c t The crack path of surface initiated rolling contact fatigue was investigated numerically based on the asperity point load mechanism. Data for the simulation was captured from a gear contact with surface initiated rolling contact fatigue. The evolvement of contact parameters was derived from an FE contact model where the gear contact had been transferred to an equivalent contact of a cylinder against a plane with an asperity. Five crack propagation criteria were evaluated with practically identical crack path predictions. It was noted that the trajectory of largest principal stress in the uncracked material could be used for the path prediction. Different load types were investigated. The simplified versions added some understanding but the full description with cylinder and asperity pressures was required for accurate results. The mode I fracture mechanism was applicable to the investigated rolling contact fatigue cracks. The simulated path agreed with the spall profile both in the entry details as in the overall shape, which suggested that the point load mechanism was valid not only for initiation but also for rolling contact fatigue crack growth. Ó 2011 Elsevier Ltd. All rights reserved.

1. Introduction Many mechanical applications contain components that interact repeatedly with pure rolling or rolling with small relative sliding. Some examples are gears, bearing, cams or even the wheel-rail contact of trains. When the rolling contact loads are high, rolling contact fatigue (RCF) may become life limiting for the structural components. RCF damage is characterized by cracks and small craters in the contact surfaces. Typical examples are presented in Fig. 1 [1]. Depending on the damage size one can distinguish between micro- and macro-scale contact fatigue. Surface distress (SD) is widely used to designate micro-scale contact fatigue. The damage has then a size comparable to the dimensions of asperities on the contacting surfaces. Macro-scale contact fatigue is commonly designated as spalling or pitting. Here the nomenclature of spalling by Tallian [2] is used. It may be noticed that a spall is the chipped off material, which by removal leaves the spalling craters. However, the crater may also be referred to as a spall. The overall gear picture in Fig. 1a shows spalls in pinion or driving teeth flanks. The spalls are located along the pitch line or roll circle with initiation points at approximately the same longitudinal position in the dedendum. Fig. 1b and c shows top and cross-sectional views of typical spalls. Talysurf measurements were performed on three spalls on the gear wheel, see Fig. 1d, where the x-axis indicated the rolling direction. In order to compare the spall profiles with numerically predicted crack paths, the curvature of the gear wheel was substracted and the spall profiles were translated to a common initiation point.

⇑ Corresponding author. Tel.: +46 8 7907667; fax: +46 8 4112418. E-mail address: [email protected] (B. Alfredsson). URL: http://www.kth.se (B. Alfredsson). 0013-7944/$ - see front matter Ó 2011 Elsevier Ltd. All rights reserved. doi:10.1016/j.engfracmech.2011.07.012

D. Hannes, B. Alfredsson / Engineering Fracture Mechanics 78 (2011) 2848–2869

2849

Nomenclature a crack length al, ap contact half-width, contact radius c crack size in y-direction E Young’s modulus hasp asperity height KI, KII mode I and II stress intensity factors KI,cl closure limit p0 l cylindrical maximum Hertzian pressure spherical maximum Hertzian pressure p0p pFE, pHertz FE and Hertzian contact pressure Pl, Pp normal line force, normal point force qFE, qHertz FE and Hertzian traction spherical maximum Hertzian tangential traction q0p Qp tangential point force r, rasp, rpl sphere radius, asperity radius, plastic zone size R radius of curvature or R load ratio Ra average roughness Rp, Rv maximum profile peak height and valley depth t thickness of FE model T T-stress xc, xd position of initial crack, position of cylindrical load x, y, z cartesian coordinates b, bentry crack angle, entry angle relative to contact surface surface distress crack angle bSD DKI mode I stress intensity factor range DKI,eff mode I effective stress intensity factor range DKth fatigue threshold value root mean square wavelength kq m Poisson’s ratio lasp asperity coefficient of friction rx, rz, sxz cartesian stresses rN, rC, rT stress normal, colinear and tangential to crack faces rR residual surface stress rY,mt, rY,mc monotonic tension and compression yield stresses rY,cycl cyclic yield stress rh hoop stress sY,cycl cyclic yield stress in simple shear  average value 0 initial value max, min maximum and minimum value

1.1. Mechanisms – previous work The literature contains numerous observations on the RCF damage process and the influence of interacting parameters. The early work by Way [3] in 1935 introduces the damage and the failure atlas by Tallian [2] displays examples of intermediate damage stages. Ding and Rieger [4] have published scanning electron micrographs of cracks and craters. The damage process consists of three stages. Firstly, inclined surface micro-cracks initiate. Secondly, some of these grow into the material, gradually turning towards a surface parallel path. The depth of this spalling bottom often corresponds to the depth with maximum effective stress. Finally, a piece of material separates from the surface giving the craters in Fig. 1. Further examination reveals two different mechanisms, which differ by initiation at the surface or below it. Both mechanisms lead to surface craters and agree in the later part of crack propagation. Surface initiated craters often display a v- or sea-shell shape, Fig. 1b, with the apex directed against the rolling direction. According to Tallian [2] the entry angle, bentry < 30°. Smaller ranges were reported by Bastias et al. [5], bentry = 20–24° and by Dahlberg and Alfredsson [1], bentry = 25–30°. The crack profiles in Fig. 1d have bentry = 23–30°. The sub-surface initiated craters are however irregular in shape and lack the shallow entry angle of the surface initiated spall. Here Tallian [2] reports bentry > 45°.

2850

D. Hannes, B. Alfredsson / Engineering Fracture Mechanics 78 (2011) 2848–2869

Fig. 1. Rolling contact fatigue damage in a gear wheel.

Many researchers have investigated lubricated rolling cylinder contacts in order to explain crack growth. However, the cylinder contact gives compressive stresses on all planes and closure of any prospective crack. The first explanation for surface initiated RCF was proposed by Way [3]. He suggested that the lubricant somehow enters an inclined crack; when overrolled a lubricant pressure introduces a mode I crack tip opening. This mechanism has over the years been refined by many researchers, but still it exhibits two major drawbacks. Firstly, it cannot explain crack initiation. Secondly, Bower [6] found through numerical simulations that lubricant pressure inside the crack can only give the damage profile through a complex load sequence. Keer and Bryant [7] suggest crack growth due to large cyclic mode II crack tip loads. The major argument against this mechanism is the experimental difficulty to produce mode II or mixed-mode I and II fatigue crack growth without crack arrest or kinking to mode I, see Bold et al. [8]. 1.2. The asperity point load mechanism Olsson [9] proposes the asperity point load mechanism, where the nominal stress field due to the two-dimensional rolling cylindrical contact is locally disturbed by the presence of three-dimensional asperity contacts. The local presence of asperities explains that RCF develops at distinct positions, see Fig. 1a. A two-dimensional mechanism can only explain line damage across the complete contact width. When a surface asperity enters the rolling cylindrical contact, the local point contact load gives tensile and forward directed surface stresses, which suffice to initiate fatigue [10]. Note that the asperity mechanism solves the classic challenge of displaying tensile stresses that are required to open a prospective fatigue crack. Fracture mechanics [11] shows that cracks from stationary but repeated point contacts follow the direction of largest mode I. The geometrical similarities between these point contact cracks and early rolling contact fatigue cracks suggest that they are the result of the same mechanism and that mode I is crack driving also at RCF. This paper describes the propagation path of surface initiated RCF cracks using the asperity point load mechanism. The purpose was to predict the crack path using measured properties from gear contacts that had suffered from surface initiated RCF. The hypothesis was that RCF crack growth followed ordinary fatigue crack laws, but was due to complex stress cycles. Correct path prediction is a necessary requirement for the asperity point load model to hold as an explanation for RCF crack growth. The work utilized earlier results [10] that predict fatigue initiation in front of the same gear tooth asperity and that model the plastic deformation with material hardening during early cycles.

D. Hannes, B. Alfredsson / Engineering Fracture Mechanics 78 (2011) 2848–2869

2851

The study of RCF damage in gears was based on two separate computational models. Firstly, a finite element contact model was used to determine the evolvement of the load cycle, including asperity contact parameters, when the asperity was overrolled. Secondly, having separated asperity from cylindrical loads, closed-form analytical contact solutions for spherical and cylindrical contacts were used to compute the stress field below the contacts, which was utilized in a numerical crack growth model for RCF. The purpose of the numerical crack growth model was to simultaneously improve the path prediction and reduce the computation time compared to an FE based crack path simulation. 2. Problem description – background RCF is complex with a large number of influencing parameters such as: load, rolling and transverse curvature, surface roughness, shape and height of individual asperities, lubrication film, additives, contaminents, slip, rolling velocity, coefficient of friction, material hardness, inclusions, micro-structure, fatigue limit, crack growth parameters, and surface treatment, which separately and interactively control the problem. Furthermore, wear and plasticity change the surface properties. The contact situation that initiated damage may not remain when the spall has developed some million cycles later. Due to the complexity, simplifications were made and the analysis was focused on some important parameters. The geometry including the surface asperity and relative movements was simplified. Wear and lubrication were only introduced through wear depth and lubrication layer. However, the simplified problem contained the important geometrical, material and load properties of the asperity point load mechanism for RCF. 2.1. Geometric model Spalling in the pinion in Fig. 1a was used as a case study. When the teeth flanks of the driving (pinion) and the faces of the driven (follower) wheel interact, the contact geometry results in a difference in velocity and slip between the contacting surfaces. On the pinion, slip increases and passes through zero at the pitch line, where pure rolling is observed. Thus, slip is negative in the zone below the pitch line where RCF initiates. The radius of curvature of the involute gear profile in Fig. 1a was assumed constant over the gear teeth and equal R = 13.6 mm [1]. The gear contact was therefore modelled as a cylinder with radius R/2 against a plane, see Fig. 2a. A small axisymmetric asperity was introduced on the plane surface, which represented the pinion. Its position corresponded to that just below the pitch line, where spalls usually initiate [2], see Figs. 1a and 2b, i.e. where slip is small and negative. The asperity had a cosine shaped cross-section with height hasp and radius rasp. The dimensions were based on 30 Talysurf measurements by Dahlberg and Alfredsson [1] on the teeth of the gear wheel in Fig. 1. Some representative profiles are displayed in [1]. The average roughness, Ra = 0.730.85 lm, the maximum profile peak height Rp = 3.13.8 lm, the maximum profile valley depth Rv = 3.26.4 lm and the root mean square wavelength, kq = 35110 lm [10]. According to the simulations by Dahlberg and Alfredsson [10] an asperity with 2 lm original height deformed during first over-rolling to 1.5 lm where it remained unchanged at subsequent over-rolling. This height of the asperity corresponded to the height after running-in, i.e. after plastic deformations flattened the asperity. The position of the cylinder centre with respect to the asperity centre was introduced as xd. The fracture mechanical investigation required an initial surface crack, which was assumed to be straight and situated at a distance xc from the centre of the asperity, see Fig. 2. The position xc was selected to correspond to the position of the maximum tensile surface stress in front of the point loaded asperity. The earlier work [10] predicts fatigue initiation at this position with the Findley criterion and material parameters from uni-axial experiments. The crack length and the angle to the contact surface were designated as a and b, respectively. A subscript naught indicated initial values. At standing contact fatigue experiments [12] on the present material the cone crack started perpendicular to the surface [11]. Since the initiation point in front of the asperity was subjected to a similar stress cycle as the ring/cone crack at standing contact fatigue, b0 = p/2

Fig. 2. Schematics (not on scale) with geometric dimensions of the asperity, spall and rolling cylinder.

2852

D. Hannes, B. Alfredsson / Engineering Fracture Mechanics 78 (2011) 2848–2869

Table 1 Summary of the introduced variables for the geometric model. rasp (lm)

hasp (lm)

R (mm)

xc (lm)

a0 (lm)

b0 (rad)

50

1.5

13.6

20

3

p 2

Table 2 Mechanical properties for the surface of the case hardened gear steel [10]. E (GPa)

m (–)

rY,mt (MPa)

rY,mc (MPa)

rY,cycl (MPa)

rR (MPa)

206

0.3

813

1803

1500

0

Fig. 3. The asperity point load model. The xz-plane is the symmetry plane.

in Table 1. The initial crack length was selected in the same range as the maximum profile valley depth, a0 = 3 lm. A short start crack length was desirable for comparison with surface distress cracks. Table 1 contains data for the geometric model in Fig. 2. 2.2. Material model During running-in of the gear, some plastic deformations were expected with isotropic hardening and kinematic shift of the yield surface. However, after running-in with elastic shake-down as was determined by Dahlberg and Alfredsson [10], isotropic linear elastic material behaviour was here assumed and motivated a posteriori. The parameters E, m and rY in Table 2 represent Young’s modulus, Poisson’s ratio and the yield stress, respectively. The monotonic and cyclic yield data are included into Table 2 as references. Gear teeth are usually case hardened resulting in graded material properties, i.e. the properties of the case differ from those of the core. The investigated gear material followed Swedish Standard 2506 and after heat treatment, the material properties were approximately constant to a depth of 0.6 mm [10]. Thus, non-graded surface material data was used. The homogenity assumption was valid as long as the cracks remained close to the surface. The case hardening also resulted in a compressive biaxial residual surface stress, rR. However, after running-in of the gear, the residual surface stress was close to zero in the rolling direction [1], hence rR was taken as 0. 2.3. Asperity point load model The gear contact in Fig. 1a was modelled as a cylindrical Hertzian pressure distribution with maximum pressure p0l = 2.39 GPa [10] and contact half-width al, see Fig. 3. Lubrication separates the gear contact surfaces. According to elastohydrodynamic (EHD) lubrication [13], the thickness of the lubrication film will normally exceed 0.1 lm. It ensures normal load transfer, but removes the tangential one. Hence, frictionless full film lubrication was assumed for the cylindrical contact. The three-dimensional asperity contact was modelled with a spherical Hertzian pressure distribution with maximum pressure p0p, tangential traction q0p and contact radius ap. The asperity was assumed to break through the lubrication film. Hence, tangential traction was transmitted over the asperity contact with the coefficient of friction lasp = 0.3. Since the slip was negative at the asperity position, the tangential asperity load was directed opposite to the rolling direction. Note that the magnitude of the spherical contact parameters depended on the position of the cylindrical contact, the contact geometry in Table 1 and the material properties in Table 2. 3. FE contact simulation The finite element (FE) program ABAQUS (6.9) was used to numerically determine the evolvement and distribution between the cylindrical and spherical contact loads during an over-rolling load cycle.

D. Hannes, B. Alfredsson / Engineering Fracture Mechanics 78 (2011) 2848–2869

2853

Fig. 4. FE geometry and mesh refinements.

3.1. FE model The FE model consisted of two three-dimensional bodies: a right cuboid solid with an asperity on its otherwise flat upper contact surface and a section of a cylindrical roller, representing the pinion and follower, respectively, see Fig. 4a. The thickness of the FE model was t = 1/3 mm. Due to symmetry, only half of the geometry was modelled and plane strain was enforced throughout the y-direction, i.e. at both y = 0 and y = t. The dimensions of the model were larger than approximately ten times a characteristic dimension of the contact zones: al in the x and z-directions and ap in the y-direction. Comparison with a three times thicker model showed less than 4% difference in the contact pressure. The mesh in Fig. 4a contained 73,949 eight-node brick elements with several mesh refinements at the contact surfaces and around the asperity. Fig. 4b represents a close-up view around the asperity. Three-dimensional transition elements were used for the 3-refinement [14]. The element size at the centre of the asperity in Fig. 4b was approximately 1.4 lm. The FE simulation consisted of two steps. Firstly, the roller was indented into the flat surface with a prescribed vertical displacement in order to get a maximum contact pressure equal to p0 l. Secondly, the roller was subjected to a prescribed rotation and the cuboic body to a horizontal translation. The rotation and the horizontal translation were related through a constraint equation, which ensured rolling with 3.7% negative slip. Large deformations kinematics were used. The contact interaction between the cylindrical roller surface (master surface) and the upper surface of the cuboic solid (slave surface) was defined as hard normal contact, thus no penetration was allowed, combined with a Coulomb friction law for tangential interaction. The coefficient pffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiof friction between the surfaces was defined through a user subroutine as lasp for nodes on the asperity, i.e. when x2 þ y2 6 ap , and zero for contact on the flat surface. The contact constraints were enforced using the augmented Lagrange method. 3.2. Hertzian contact parameters The total FE contact pressure pFE and traction qFE at y = 0 were approximated with cylindrical and spherical Hertzian contact pressures at each time step of the FE load cycle. The total contact pressure pFE, was extracted along two separate positions. Firstly, along y = t, where the contact pressure was assumed independent of the asperity contact and allowed computation of p0 l. Secondly, pFE was extracted along y = 0, where the asperity contact pressure profile was isolated by substraction of the cylindrical contribution and then p0p was computed. Furthermore, the FE contact traction was translated to q0p. With gross sliding on the asperity, the magnitude of q0p was recovered by lasp pHertz(x = 0). Note that when the asperity was situated in the cylindrical contact zone, both p0p and the cylindrical contact pressure at x = 0 contributed to pHertz(x = 0). The contact loads were computed for an uncracked substrate, as the increased compliance due to a propagated RCF crack was assumed to have a negligible effect on the contact loads. It was noted that a sphere with radius r  4 mm would meet the long crack in Fig. 1d at x = 1.5 mm and z = 0.3 mm. A comparison with the corresponding Hertzian contact with ap/r  0.02/ 4 = 0.005 verified the assumption. The total normal pressure pHertz was obtained for different values of xd, i.e. for different positions of the rolling contact, by superposition of the cylindrical and spherical Hertzian pressures. In Fig. 5a contact between the asperity and the cylindrical surface was established, but the asperity was still situated outside the cylindrical contact zone. The Hertzian formulation of the asperity contact was assumed to be centered on the geometrical asperity, which lead to a small offset between the FE and

2854

D. Hannes, B. Alfredsson / Engineering Fracture Mechanics 78 (2011) 2848–2869

Hertzian asperity pressures. In the FE model the contact occured first on the flank of the asperity and not precisely at its centre. Fig. 5b corresponds to the critical position, where the asperity entered the cylindrical contact zone and where the surface stress rx in front of the asperity became maximal for the Hertzian stress formulation. The difference in maximum pressure between the two models was approximately 7% in Fig. 5b. The geometry of the asperity caused a reduction of pFE just outside the asperity contact zone, which pHertz could not account for, see Fig. 5c. Fig. 5d corresponds to the cylindrical contact centered upon the asperity. Hertzian theory [15] gave the dimensions of the cylinder and asperity contacts,

2Rð1  m2 Þ p0l ; E ! 2 pr2asp 2Rð1  m Þ ; p0p ap ¼ E 4r 2asp þ p2 Rhasp

al ¼

ð1aÞ ð1bÞ

where the relative curvature between the contacting surfaces was considered. The mean curvature at the centre of an axisymmetric cosine shaped asperity was used for the asperity contact. During the RCF load cycle ap varied, as opposed to al, which remained constant. For p0l = 2.39 GPa, al = 287 lm. The normal and tangential asperity contact loads were assumed to have the same contact radius. The evolvement of p0p and q0p during a load cycle are given in Fig. 6, where a moving average filter was applied to reduce variations between consecutive positions of the cylindrical contact. During over-rolling p0p and q0p remained fairly constant, which through Eq. (1b) yielded a rather constant ap 6 22.5 lm. For the case study of Fig. 1a, the extreme values of the maximum pressure and traction were 5.0 GPa and 2.1 GPa, respectively. Outside the interval 363.2 lm < xd < 363.6 lm the asperity load was zero and the fatigue crack was only subjected to the cylindrical contact load.

Fig. 5. Hertzian and FE pressure p and traction q at y = 0, when the asperity progessively entered into full contact (xd = 0 lm).

D. Hannes, B. Alfredsson / Engineering Fracture Mechanics 78 (2011) 2848–2869

2855

Fig. 6. Evolvement of the maximum asperity pressure and traction. A moving average filter was applied.

3.3. RCF initiation region Fig. 7a and b shows pFE and the surface stress rx in the vicinity of the asperity for the cylindrical contact position that yielded max(rx(x, xd)), i.e. the maximum FE surface stress rx, for all xd and positions x along the symmetry line (y = 0). The results show that max(rx(x, xd)) occured near the border of the asperity contact. Observe also that a large zone in front of the asperity was not subjected to any contact pressure and experienced tensile surface stress rx. Fig. 7a illustrates that the asperity pressure was spherical with a circular radius ap.

(a)

(b) Fig. 7. (a) FE results for contact pressure pFE and (b) surface stress rx in the vicinity of the asperity for xd = 260.6 lm. The white semi-circle represents the asperity size and the position max(rx(x, xd)) is indicated with .

2856

D. Hannes, B. Alfredsson / Engineering Fracture Mechanics 78 (2011) 2848–2869

Fig. 8. Comparison of FE results and closed-form stress solutions based on Hertzian contacts with LT 23. The position of max(rx(x, xd)) is indicated with

.

The FE surface stress rx(x, xd) was extracted along the symmetry line. Fig. 8a presents maxðrx ðxÞÞjxd , i.e. the maximum surface stress rx for varying x and fixed xd. Note that the tensile surface stresses induced by the asperity load were the highest when the asperity entered the cylindrical contact. Hence, max(rx(x, xd)) = 2.1 GPa and occured for xd = 260.6 lm, see Fig. 7b, and was located at x = 28.1 lm < rasp. Fig. 8b shows max(rx(xd))jx, i.e. the maximum surface stress rx for varying xd and fixed x. The FE geometry was not an infinite body, which explained non-zero rx when the asperity was not loaded by the cylinder. The zone with high tensile surface stress in front of the asperity, i.e. when x > 0, was identified as a potential RCF initiation region. An initial crack perpendicular to the surface would be subjected to rx. 4. Numerical crack growth model for RCF The RCF crack path was modelled in the symmetry plane (y = 0). The rolling contact with an asperity is three-dimensional, but by considering the symmetry plane through the asperity, a two-dimensional crack view was motivated and two-dimensional fatigue crack growth was investigated for the symmetry plane, see Fig. 3. This simplification is valid for spur gears and fully lubricated helical gears if the transverse friction component for the asperity is neglected. Note that by transforming the FE contact pressures to Hertzian pressures, closed-form expressions exist for the stresses at all positions in the assumed half-plane [15,16]. The stress state from the asperity point load was superimposed onto the twodimensional stress state from the cylindrical contact. The resulting stress field was used to model crack propagation with Matlab (R2009b). 4.1. Load models At this point, the gear contact with an asperity on the pinion surface has been transformed into a moving cylindrical contact pressure, fixed spherical contact pressure and fixed spherical contact traction both with varying magnitudes, see Fig. 6. Alternatively, the pressures could be integrated into line and point forces. The cylindrical contact was represented by a concentrated normal line force

Pl ¼

pp0l al 2

:

ð2Þ

Similarly, the asperity contact was characterized with a concentrated normal point force

Pp ¼

2pp0p a2p ; 3

ð3Þ

and a tangential point force

Qp ¼

2pq0p a2p ; 3

which was directed opposite to the rolling direction, see Fig. 3, as slip was negative and thus q0p < 0.

ð4Þ

D. Hannes, B. Alfredsson / Engineering Fracture Mechanics 78 (2011) 2848–2869

2857

Table 3 Definition of the six loading types (LT). Load case

Concentrated loads

Distributed loads

Pp Pp and Qp Pp, Qp and Pl

LT 11 LT 12 LT 13

LT 21 LT 22 LT 23

The closed-form expressions for the stress fields are available as elementary contributions [15,16]. For completeness, these are presented in Appendix A. Three combinations of elementary loads were investigated. Firstly, the load consisted only of the normal asperity load (Pp). Secondly, the tangential component on the asperity was included (Pp and Qp). Thirdly, both asperity loads were superimposed onto the normal cylindrical load (Pp, Qp and Pl). For all three load cases, both concentrated forces and distributed pressures were considered, see Table 3. Fig. 8 shows the surface stress rx along the symmetry line, based on the elementary solutions for LT 23. Its absolute maximum was equal to 3.2 GPa, which was larger than the FE maximum at 2.1 GPa, and occured for xd = 266.6 lm, see Fig. 8a. Its position on the symmetry line was x = 20.3 lm, see Fig. 8b and motivated xc in Table 1. For the numerical solutions, x and xd were discretized, with increment sizes comparable to those used in the crack growth model. The jagged appearance of the closed-form solution based on Hertzian contacts in Fig. 8a is due to the dicretization in x and xd. The surface stress rx based on closed-form expressions in Fig. 8a agreed qualitatively with the finite element results. However, differences in magnitude and position were observed. The closed-form solution contained a higher max(rx(x, xd)). It occured earlier in the roll cycle and closer to the asperity centre. The closed-form solution in Eq. (A.4a) has the maximum at the contact boundary and decays with 1/x2 for increasing x. It gives an upper bound for the maximum surface stress. The FE solution suffers from the mesh discretization that gives the maximum value at least one element outside the boundary. With ap = 22.5 lm and 1.4 lm element side the stress maximum will decrease about 12% from this effect. Furthermore, the eight-node brick elements use linear shape functions that are less optimal in resolving steep stress gradients. For both solutions, the maximum surface stress rx occured just outside the cylindrical contact zone, i.e. at approximately the position x = xd + al, which illustrated the influence of the position of the contact boundaries. 4.2. Crack growth model During the non-proportional fatigue cycle the crack was subjected to both normal load (mode I) and in-plane tangential load (mode II). Hence, plane mixed-mode crack propagation was expected in the symmetry plane with a curved crack path. The crack path was predicted assuming that the mode I crack growth mechanism was dominant. Moreover linear elastic fracture mechanics (LEFM) was assumed applicable. The stress intensity factors (SIF) KI and KII were calculated with the superposition principle, as the load on the external boundaries of a structure with an unloaded crack can be translated to an equivalent distributed stress field on the crack faces. The stress normal rN and tangential rT to the crack path in the uncracked half-plane gave the load on the crack faces. The SIFs for a given crack length a were computed by integration of the elementary formulae for line forces normal and tangential to the crack faces of a through-thickness straight edge crack [17]. The SIFs were expressed by weighted integrations along the crack faces of the normal stress,

K I ðaÞ ¼ 2

rffiffiffiffi Z a

p

1

0

 5 rN ðnÞ  pffiffiffiffiffiffiffiffiffiffiffiffiffi ffi 1:3  0:3nð4Þ dn; 2 1n

ð5aÞ

and the tangential stress,

rffiffiffiffi Z a

K II ðaÞ ¼ 2

p

1 0

 5 rT ðnÞ  pffiffiffiffiffiffiffiffiffiffiffiffiffiffi 1:3  0:3nð4Þ dn: 2 1n

ð5bÞ

The integrands in Eq. (5) were singular at the crack tip, i.e. for n = 1. To ensure accuracy of the numerical integration the density of integration points near the crack tip was increased using a non-uniform distribution based on a double sine function. At least 500 integration points were used along the crack faces with linearly interpolated stresses. The ratio between the maximum and minimum spacing between integration points was approximately 1.6  108. Various crack path prediction criteria (CPPC) exist in the literature [18,19]. Five such mode I criteria were implemented, see Table 4. The criteria give a prediction of the deflection angle, which corresponds to the change in the crack angle b for the next crack increment. The sign of the deflection angle for the SIF based criteria is given by the sign of KII. The load cycle was non-proportional in mode I and II. Therefore a criterion was required to select the instance during the load cycle when the crack path direction should be determined. Firstly, it was noted that KI was negative for the major part of the load cycle and that a closure value, KI,cl, was required. For simplicity KI,cl = 0 was assumed. Secondly, the fracture surface in Fig. 1b and surface profiles in Fig. 1d were investigated. The surfaces show an irregular profile without wear marks, which suggested that limited crack face sliding had occurred during crack growth. Hence, it was assumed that when KI < KI,cl = 0 the crack faces were fixed, KII did not change and the crack did not grow. Thirdly, the remaining part of the load cycle contained the critical

2858

D. Hannes, B. Alfredsson / Engineering Fracture Mechanics 78 (2011) 2848–2869 Table 4 Description of the five crack path prediction criteria (CPPC). Criterion

Description

CPPC CPPC CPPC CPPC CPPC

Maximum tangential principal stress criterion Maximum tangential stress criterion Minimum strain energy density criterion by Sih Maximum energy release rate criterion by Nuismer Empirical criterion by Richard

1 2 3 4 5

instance that controlled the crack growth direction. The instance with the maximum KI was assumed to control the crack growth direction. An alternative critical instance occured when the hoop stress rh at the crack tip was the maximum given that KI > KI,cl = 0. Some crack path prediction criteria assumed positive KI, which was verified a posteriori. The maximum tangential principal stress criterion assumes crack propagation perpendicular to the largest principal stress direction in the uncracked material position that corresponds to the crack tip. This criterion can also be called the zero shear stress criterion. The maximum tangential stress criterion proposed by Erdagon and Sih [20] suggests a radial crack extension from the crack tip in the plane perpendicular to the direction of largest hoop stress. This criterion is also known as the maximum hoop stress criterion for which an explicit expression of the deflection angle can be formulated in terms of the SIFs, when dominance of the first order stress field at the crack tip is assumed. Another criterion is based on the work by Sih [21]; the crack propagates in the direction where the strain energy density presents a stationary minimum. The minimum strain energy density criterion by Sih is also called the S-criterion, in order to dissociate it from other strain energy density based criteria. A characteristic of these criteria is the path dependence on Poisson’s ratio. Palaniswamy and Knauss [22] proposed a criterion based on the energy release rate, which predicts crack propagation in the direction corresponding to the maximum value of the energy release rate. Different plane mixed-mode energy release rates in a homogeneous isotropic elastic material have been proposed, among which the formulation by Nuismer [23]. The maximum energy release rate criterion by Nuismer gives the total energy release rate in terms of the SIFs at the crack tip of the propagated branch, which were determined by using a continuity assumption for the stress field when the crack extension of a deflected crack tends to zero. Finally, the criterion by Richard [24] is an empirical criterion that expresses the deflection angle as function of the SIFs at the crack tip. 4.3. Calculation model – method The numerical implementation of the crack growth model was based on a double discretization. Firstly, the Hertzian contact pressures and traction were defined for a series of discrete positions xd of the cylindrical load. These xd ranged from 500 lm to 450 lm, with an increment equal to 2.85 lm. In order to increase accuracy of the numerical crack growth model, this increment was reduced to 0.285 lm with linear interpolation of the data. Also the load cycle was extended to xd = 2 mm with load steps that only contained the cylindrical pressure using larger increments. Secondly, the crack path was discretized with a minimum crack increment of 0.1 lm. The crack increment was increased for a > 0.1 mm, in order to decrease the total computational time. It was assumed that during each load cycle the crack would propagate with KI,max > KI,cl = 0, which was enabled by setting the threshold value for fatigue crack growth DKth to zero. 5. Crack growth results for RCF Firstly, the influence of loading type was investigated, which added understanding of the asperity point load mechanism for RCF. Secondly, the five CPPC were compared with the objective to select the one most suited for the continued investigation. Finally the evolvement of the SIFs was studied. 5.1. Influence of loading type on crack path The influence of loading type was studied using the maximum tangential principal stress criterion. Fig. 9a shows the predicted crack paths superimposed on an experimental spall profile from Fig. 1d. The x and z-axes use equal scales. The experimental spall profiles were translated 5 lm vertically to take mild wear around the pitch line into account [25]. Normal asperity load alone (LT 11 and 21) resulted in a crack steeper than the experimental spall profile; adding the tangential asperity load (LT 12 and 22) resulted in an even steeper crack path prediction. Thus, the asperity load alone predicted to deep overall fatigue cracks in Fig. 9a. When the normal cylindrical load was added, the predicted crack path became more shallow and close to the spall profile. When the paths from each pair of concentrated and distributed formulations were compared, the one from the distributed contact load always predicted the steeper crack path. In fact, the distributed load that included normal cylindrical load as well as normal and tangential asperity load (LT 23) was consistent with the measured spall profile. The initial crack path predictions in Fig. 9b displayed slightly different behaviours. Fig. 9c shows KI,max during a load cycle at each crack length. After an initial increase for short crack lengths, KI,max decreased rapidly for increasing crack lengths to relatively low, but positive KI,max. The tangential asperity load resulted in

D. Hannes, B. Alfredsson / Engineering Fracture Mechanics 78 (2011) 2848–2869

(a)

(b)

(c)

(d)

2859

Fig. 9. Influence of LT with CPPC 1 on (a) crack path, (b) early crack path, (c) KI and (d) KII when KI is the maximum.

increased KI,max, compared to the normal asperity load alone. The opposite effect, but to a lesser extent, was obtained when the cylindrical load was added. Finally, distributed contact loads always predicted lower KI,max than concentrated loads. The KII value when KI was the maximum, i.e. K II jK I;max was a magnitude lower, see Fig. 9d. After an initial increase for short crack lengths, also the magnitude of K II jK I;max decreased rapidly for increasing crack lengths. 5.2. Influence of loading type on surface stresses Fig. 10 displays the crack opening stress, rN, and in-plane shearing stress, rT, at two crack positions and the over-rolling cycle was represented by xd. By investigating these stresses, some insight was gained into initiation properties of the point load mechanism for RCF. One crack position was (x, z) = (xc, 0), which was the expected crack initiation position at the surface. The other was the initial tip position at (xc, a0) for the crack path simulations. At both these positions, rN = rx and rT = sxz, since the initial crack tip angle was b0 = 90. In Fig. 10a and b, the stresses are non-zero when the asperity is in contact or when the position xc is inside the asperity contact zone. When the mid-point of the cylinder was at smaller or larger xd, both the asperity and the crack position were unloaded by the cylindrical contact, which still contributed to the stress state under the surface. The importance of asperity pressure and traction as well as cylindrical pressure was illustrated by LT 11–23. In Fig. 10a, LT 11 illustrates how the point force on the asperity creates a tensile rN on the surface. When the normal asperity load is refined to the asperity pressure with LT 21, rN decreases substantially to compressive values when xc is inside the asperity pressure, i.e. ap > xc. However, there remain two tensile stress peaks at the beginning and end of the load cycle. Adding the friction on the asperity with traction increased rN substantially for all xd, but did not change the principal shape of the stress profile. The peak rN occurred when ap = xc. The two last curves illustrate the influence of the two-dimensional load as a line force, LT 13, or a cylindrical pressure, LT 23. The results for LT 13 were excluded for 266.6 lm < xd < 266.6 lm because of the stress singularity. The line force description did not add any information for rN in Fig. 10a and the cylindrical pressure description

2860

D. Hannes, B. Alfredsson / Engineering Fracture Mechanics 78 (2011) 2848–2869

Fig. 10. Stress results near the surface obtained with the different loading types (LT).

again lowered rN at over-rolling to compressive values and decreased the second stress peak. The small double peak at contact exiting was the result of xc exiting first the asperity pressure and then the cylindrical pressure. Finally, each over-rolling gave two stress cycles with approximately the same stress amplitude and peak stresses. This effect diminshed rapidly below the surface, as the second tensile stress peak disappeared at z = 4.6 lm. Therefore, the second peak was not included in the crack growth simulations. The total stress state also included rT – 0 on the surface from friction on the asperity; see the cycle for LT 22 and 23 in Fig. 10b. At 3 lm below the surface, the normal stresses in Fig. 10c have decreased substantially. In particular the peak for LT 23 had almost disappeared. The shear stress in Fig. 10d, on the other hand, has increased substantially for all loading types. Thus, the principal stresses rapidly changed both magnitude and direction below the surface. The substantial shear stress suggests that already for very short cracks the largest mode I direction was not longer parallel to the contact surface. In particular the shear stress value when rN was large in Fig. 10a and c was expected to influence the direction of the first increment. The crack growth prediction was therefore expected to start with a rapid turn. Note that LT 12 and 13 display positive rT(xc, a0) and should therefore start with turning backwards, which they do in Fig. 9b. The other paths directly turned in the forward or rolling direction following the negative rT value in Fig. 10d.

5.3. Influence of crack path prediction criterion The five CPPC in Table 4 were investigated with LT 23, see Fig. 11. The initial crack length was a0 = 0.1 lm to capture any potential early crack kink prediction by any of the CPPC. The different CPPC gave close to identical crack paths, see Fig. 11a. The difference in crack depth, i.e. z, at a given crack length, between the crack path of the maximum tangential principal stress criterion (CPPC 1) and the other criteria was less than 1%, see Fig. 11b. Close to identical crack paths yielded close

D. Hannes, B. Alfredsson / Engineering Fracture Mechanics 78 (2011) 2848–2869

2861

pffiffiffiffiffi to identical KI,max, see Fig. 11c. The maximum KI,max was equal to 3.45 MPa m and occured for a crack length of 1.6 lm. The K II jK I;max was significantly smaller and differed between the criteria, see Fig. 11d. Note that the results in Fig. 11d were trunpffiffiffiffiffi pffiffiffiffiffi cated. Indeed the maximum tangential principal stress criterion (CPPC 1), varied in the range 36 kPa m to 11 kPa m, pffiffiffiffiffi pffiffiffiffiffi whereas the other criteria oscillated between 22 kPa m to 2 kPa m. For the latter criteria numerical stability problems were observed. The maximum tangential principal stress criterion, on the other hand displayed stable crack path predictions throughout the simulation.

5.4. The asperity point load model Focusing on LT 23 with CPPC 1, the crack path prediction with the critical growth instance defined by KI,max was compared to the profiles of three spalls in Fig. 12a. Moreover the crack path predicted with the propagation instance determined by the maximum hoop stress at the crack tip occuring when no closure was present, i.e. at rh;max ¼ maxðrh ðh; xd ÞjK I >K I;cl ¼0 ), was superimposed. Crack growth at rh,max did hardly modify the predicted crack path: the largest relative difference in depth between the two predicted crack paths in Fig. 12a occured at a = 2 mm and was approximately equal to 1.5%. The critical xd position corresponding to rh,max oscillated around the more stable one that gave KI,max. Fig. 12b shows the in-plane stresses at the crack tip for KI,max, where rC is colinear to the crack. Starting at 3.1 GPa, rN decreased rapidly for increasing crack lengths. At a = 10 lm, rN = 0 and remained slightly negative for increasing crack lengths. For a < 2 mm, rT < 0. At a0, rT = 300 MPa, but when the crack propagated according to the maximum tangential principal stress criterion, rT  0, which was consistent with the propagation rule. The stress rC < 0 for a < 2 mm with the minimum 607 MPa at a = 3.1 lm. Fig. 13a and b presents the non-proportional KI and KII in terms of a and xd. The gray contour levels correspond to a zero SIF. The evolvement of the KI and KII during the roll cycle was extracted for 6 selected crack lengths, see Fig. 13c and d. The x

(a)

(b)

(c)

(d)

Fig. 11. Results obtained with LT 23 and the different crack path prediction criteria (CPPC). (a) Crack path and (b) depth position of CPPC 2–5 relative to z of CPPC 1. (c) Maximum mode I SIF and (d) corresponding KII as function of crack length.

2862

D. Hannes, B. Alfredsson / Engineering Fracture Mechanics 78 (2011) 2848–2869

(a)

(b)

Fig. 12. Results for LT 23 with CPPC 1. (a) The crack path prediction using different critical instances superimposed to three spall profiles from Fig. 1d. (b) Stresses at the crack tip for KI,max.

(a)

(b)

(c)

(d)

pffiffiffiffiffi Fig. 13. SIF results for LT 23 with CPPC 1. (a) KI and (b) KII during over-rolling and crack growth, in MPa m, see text for details. Evolvement of (c) KI and (d) KII during the roll cycle for six different crack lengths: a = a0, 50 lm, 0.2 mm, 0.5 mm, 1 mm and 1.5 mm.

D. Hannes, B. Alfredsson / Engineering Fracture Mechanics 78 (2011) 2848–2869

2863

position of the crack tip as function of a is represented as a white dashed line in Fig. 13a and b and as white squares in Fig. 13c and d. Fig. 13a shows that KI was positive only in two regions, one where the asperity entered and a second where pffiffiffiffiffi pffiffiffiffiffi the asperity exited the contact. In these regions, max(KI) equalled 3.4 MPa m and 2.5 MPa m, respectively. When the asperity entered the contact zone, it yielded a positive mode I SIF for a < 2 mm. When exiting the contact zone, KI > 0 only pffiffiffiffiffi for a < 15.2 lm. In the region with a negative KI, a zone with KI lower than 80 MPa m could be observed. KI,min(a) occured when the cylindrical load was above the crack tip, i.e. later in the load cycle when a increases, see Fig. 13c. On the other hand KI,max always occured when the asperity entered the contact zone. Fig. 13b shows the evolvement of KII(a, xd). When a increased, the extreme values of KII occured later in the load cycle, i.e. for larger xd. KII,max occured just before the cylindrical contact passed over the tip and KII,min just after the passage, see Fig. 13d. Thus, the crack tip experienced a load sequence with KI,max, KII,max, KI,min and KII,min. It should be noted that the effective SIF ranges would be limited by closure and crack face friction. Crack closure effects were not taken into account in Fig. 13. 6. Discussion The analysis included some assumptions. Some were more straightforward than others. The asperity size with radius rasp and height hasp was related to surface roughness measurements [1] and the assumed axisymmetric cosine shape had some foundation in the Talysurf measurements. The real asperity shape may however differ in details from the cosine shape and therefore the real asperity pressure may differ in details from the spherical Hertz profile. Nonetheless, the major hypothesis remains, local asperities break through the lubrication film and create point loads with a local convex curvature that gives a tensile surface stress rx at contact entry. The initial crack was short, perpendicular to the surface and positioned at the maximum surface stress in front of the asperity. As long as the initial crack remained within the region with high surface stress in Fig. 7b, the simulated path was relatively unsensitive to xc. Other positions in this region would result in slightly shallower paths. Given the short size of the crack, the initial crack length and direction have low influence, while the crack will immediately kink to a profile quasi-identical to the one in Fig. 12a. The important result from the FE simulation was the distribution between the cylindrical and the asperity pressures including the sequence for them during the over-rolling. The continued analysis utilized that the pressure distribution at asperity over-rolling could be well approximated with a cylindrical Hertz pressure with superposition of a spherical Hertz pressure for the asperity contact, see Figs. 5 and 7a. Since closed-form expressions exist for the stresses below these pressures, the stresses were directly available for all points along any path below the surface. Therefore, the continued analysis could be refined more than 10 times compared with the densest part of the FE mesh while the computations became much faster than an FE analysis. The elementary formulae for KI and KII at straight edge cracks were used to approximate those with curved crack paths. It included two approximations with respect to crack curvature. The possible errors of these assumptions were estimated in Appendix B with respect to crack path curvature and in Appendix C for crack front curvature. Note that by transforming the stresses along the crack path to the local crack face orientation, the effect of crack path curvature was limited to overestimating KI with 35% and underestimating KII with 9%, for the circular-arc crack path in Appendix B. The crack front curvature could at the most give an equally large error, see Appendix C. However, the investigated spall front in Fig. 1b suggested that the crack front actually was rather straight. Hence, a/c was relatively close to 0 for the spall and the effect of crack front curvature on KI and KII should be well below the maximum value 41%, which appeared at a/c = 1 in Appendix C. LEFM was assumed with in-plane crack tip load described by KI and KII. The limits for LEFM are typically expressed as the crack being micro-structurally or mechanically short [26]. The material was pearlitic with martensite and some retained austenite. Pearlite lamellas are typically a few micrometer thick and in this case the lamellas were broken by martensite structures. The term mechanically short depends on the size of the crack tip plastic zone rpl relative to the crack length. Riemelmoser and Pippan [27] find that LEFM starts to break down for monotonic load on short cracks for a = 25100 rpl, where 100 is a strict limit and 25 may be more appropriate in applications. Suresh [28] gives an estimate of the monotonic plastic zone at plane strain that was adopted to cyclic conditions following Rice [29] by replacing KI,max with DKI and rY with 2rY,cycl. Hence,

a > 25r pl ¼

  25 DK I 2 : 3p 2rY;cycl

ð6Þ

pffiffiffiffiffi Here, DK I ¼ 1  5 MPa m, rY,cycl  1500 MPa from Table 2 and a > 0.3–7.4 lm. This limits on crack length agree with the overall transition length for short cracks in high strength steel, a > 1–10 lm that Suresh [28] found from the Kitagawa– Takahashi diagram. Based on mechanical and micro-structural comparisons, LEFM may be acceptable also at the initial crack length in the simulations. A more practical view is to use LEFM with KI and KII for the analysis and compare the predicted path and spalling crack shape. Since the predicted paths agree with the spall profiles in Figs. 11a and 12a it was concluded that LEFM and the mode I assumption was acceptable for the crack path simulations. The five evaluated CPPC suggested quasi-identical crack paths, although they were based on different theories and implementations. In fact it was hard to determine differences between the criteria. From the practical path simulation perspective in Fig. 11a and b, they were identical. In Fig. 11d small differences were noted, primarily between the maximum tangential

2864

D. Hannes, B. Alfredsson / Engineering Fracture Mechanics 78 (2011) 2848–2869

principal stress criterion, CPPC 1, and the other criteria. The KII differences were however judged as negligible but the oscillations that sometimes appeared for the four SIF based criteria were a clear drawback for these. Thus, the maximum tangential principal stress criterion was preferred. It gave accurate and stable solution while at the same time being straightforward to evaluate. The FE simulation with data from the gear contact in Fig. 1a resulted in a large tensile surface stress in front of the asperity, see Fig. 7b. The maximum von Mises stress in front of the asperity was approximately 3 GPa, which is higher than the monotonic yield stress in Table 2. However, the analyses focused on the steady-state solution after rapid isotropic hardening [10] to rY,cycl  1500 MPa in Table 2 and a kinematic shift in the stress space. Hence, the elastic stress range should be compared to the diameter of the cyclic yield surface 3 GPa. The high stress was limited to a small region, see rN in Fig. 12b. Initial yielding in an elastic–plastic model would be limited to the absolute surface region and was not expected to affect the substrate compliance or contact pressures. The high surface stress was sufficient for initiating fatigue in front of the asperity, which is examined in detail by Dahlberg and Alfredsson [10] for the elastic–plastic material description. The peak asperity pressure is embedded in the surface layer of the cylindrical contact with p0 l = 2.39 GPa. At plane strain the hydrostatic pressure in the surface of the cylindrical contact becomes 2(1 + m)p0 l/3  2 GPa, which reduces the risk for plastic deformation at the surface. The maximum spherical pressure on the asperity equals approximately 5.0 GPa in Fig. 6 and the total pressure on the asperity reaches 7.5 GPa in Fig. 5d. The elastic model represented steady-state conditions, which required that this load remained within the elastic shake-down limit. The situation experienced by the over-rolled asperity resembled repeated spherical indentation inside a hydrostatic pressure of about 2 GPa. Since hydrostatic pressure does not influence the risk for yielding, the remaining pressure of 5.5 GPa was evaluated with respect to shake-down at re ¼ 2p0p =3 ¼ 3rY , Johnson [15] finds that reversed peated spherical contact. For a fully plastic spherical contact that sustains p plastic deformation is only possible due to the Bauschinger effect. At repeated contact, kinematic hardening is avoided when the effective stress remains below rY,cycl. Hence, the pressure that remains after removing the hydrostatic component of the cylindrical pressure should be less than 9/2rY,cycl = 6.75 GPa for elastic shake-down below the contact. Since the remaining pressure was 5.5 GPa, the elastic shake-down assumption held, which is confirmed by the numerical results in [10]. The present analysis shows that when the crack has initiated in the surface due to the high tensile surface stress in front of the contact loaded asperity, the fatigue crack path follows ordinary fatigue crack growth laws. The crack turned continuously with the direction that gave the maximum mode I SIF during the load cycle with K II jK I;max  0, see Fig. 11d. Fig. 9b shows how the crack rapidly turns without a crack kink to the entry angle bentry = 20–30° that numerous researchers have documented [1,2,5,9]. The simulated path displays a stable entry direction, for 11 lm < a < 30 lm. LT 23 predicts bentry = 24.1°, which falls in the middle of the measured values. The overall path prediction for LT 23 is visible in Fig. 9a. The spalling crack behaviour and its agreement in entry detail and overall shape with the simulated mode I path suggests that the asperity point load mechanism is important for explaining the spalling crack path. For all crack lengths the KI,max always appeared at the same xd, see Fig. 13a; the position when the asperity entered the cylindrical contact. Fig. 13a and b shows the mixed-mode sequence of the SIFs. The sequence at each crack length a is identified by vertical lines in the figures resulting in Fig. 13c and d. Before the cylinder reaches the asperity both KI and KII are close to zero. The SIF cycle starts with KI,max after which KI decreases rapidly while KII increases. The sign of KII agrees with rT in Fig. 12b. Hence, KII > 0 would try to kink the crack down into the material, but KI < 0 counteracts this kink. As KI < 0,KII starts to decrease. Then KI increases and when KII,min < 0 is reached KI is only slightly negative. KII < 0 tries to kink the crack up towards the surface. As a increases, KII,min decreases in Fig. 13b while KI remains slightly negative. Hence, a highly negative KII would subdue a small crack closing KI and the crack should eventually kink to the surface for long a. The mixedmode crack load agrees with the irregular spalling crack paths in Fig. 12a and Fig. 1d and the presence of wing cracks in Fig. 1c. The presence of a negative T-stress is a stabilizing and kink resisting factor that halts the wing cracks. Fett [30] derives the T-stress for an edge crack in an infinite sheet: T = rN(a), using a Green’s function similar to those leading to Eq. (5). The present superposition of crack solutions suggested that the T-stress could be approximated as

T ¼ rC ðaÞ  rN ðaÞ:

ð7Þ

Fig. 12b gives rC(a) and rN(a) for the instance with KI,max. The stress rN is in fact slightly negative for a > 10 lm but since rC < rN, T < 0 and crack stabilizing. After the initial turn at 510 lm the simulated paths in Fig. 9b include a short range with stable b that differs from the  and LT 11–23. The normal point load overall path angle. Table 5 displays the average crack angle for 11 lm < a < 30 lm; b, alone predicts early crack paths at 22–26°, whereas the combined normal and tangential point loads predicts deeper paths at 36–57°. In applications, the early contact surface damage, surface distress [2], develops relatively early during the contact surface life and displays short cracks, length 1030 lm, that are inclined in the forward direction at and below the pitch line of pinions. Two groups of inclination angles to the surface are documented for the pinion [9]. At the pitch line where friction is absent bSD = 18–28°; before the pitch line with negative friction bSD = 41–50°. If asperities are worn during running-in, then the asperities and point loads should be relatively dominating during early cycles with surface distress formation. The line load should have relatively limited influence on surface distress and the agreement between the predicted angles with and without tangential traction in Table 5 and surface distress in gears was reasonable. pffiffiffiffiffi In Fig. 11c KI,max > 0 for all crack lengths. The value is however small, K I;max < DK th  5 MPa m [11]. Yu et al. [31] investigate crack growth at R = 1 in aluminium. They find decreasing DKth and increasing crack growth rate with increasing compressive minimum stress. Other researchers express crack growth with DKI,eff = KI,max  KI,cl, where KI,cl can be negative for

D. Hannes, B. Alfredsson / Engineering Fracture Mechanics 78 (2011) 2848–2869

2865

Table 5  average crack angles for 11 lm < a < 30 lm, see Fig. 9b. Simulated b, Load case

Concentrated loads (°)

Distributed loads (°)

Pp Pp and Qp Pp, Qp and Pl

22.2 36.4 21.3

26.1 56.9 24.1

R < 0. Silva [32] detects KI,cl < 0 at R = 1 for a relatively high strength steel; Romeiro et al. [33] find the same for a mild steel at 3 6 R 6 1; Pommier and Bompard [34] verify the finding for a mild steel with numerical simulations that predict negative KI,cl at 21% of KI,max when R = 1. The present crack tip load cycles resembled the described situations. The cycles in Fig. 13a include a small positive KI,max and a highly compressive minimum stress, i.e. large negative KI,min and thus large negative R values. Based on the findings in the literature for similar materials, it was concluded that crack growth was possible for the present load cycles. Dahlberg and Alfredsson [10] used the asperity point load mechanism to explain fatigue initiation in front of overrolled asperities. In the now described work, the forward directed crack path of RCF was simulated from the fatigue initiation site in front of the over-rolled asperity with LEFM and mode I crack growth. RCF does include some further properties that should be mentioned. The damage typically develops in the forward rolling direction and almost exclusively for negative slip. Hence, the slip direction plays an important role in the damage process. The present work showed that if the asperity is located below the pitch line of the pinion where slip is negative, then friction will increase the tensile surface stress component in front of the contact entering asperity. Crack initiation will be promoted by the point load stress and friction. The crack will grow and turn in the forward direction driven by the asperity point load in a similar way as the cone crack reported in [11] for standing contact fatigue. When the asperity exits this rolling contact, the point loaded asperity will give a tensile radial surface stress behind the asperity. Here friction from slip will decrease the surface stress from the point load and a backward directed fatigue crack is therefore less likely at contact exit with negative slip. If rolling continues but the slip direction is reversed as happens for the pinion gear tooth when the contact has passed the pitch line, then spalling damage in applications becomes scarce. The cracks that do develop may grow in the backward direction or being forward directed with an even more shallow entry angle. Note that the present Hertzian description is symmetric; if the slip direction is reversed, then backward crack growth would be predicted at asperity exit. However, two-dimensional rolling contacts in applications are always lubricated to avoid wear and other damage modes. At highly loaded rolling contacts the EHD contact profile develops [13] with a pressure spike at the exit of the lubricated cylindrical contact. This profile is numerically hard to model and was outside the present scope while it does not alter the entry contact profile of the cylindrical contact pressure. The asymmetry effect of the EHD pressure spike can however be qualitatively envisioned. The spike enhances the line load at contact exit while the pressure on the lubrication film breaking asperity is unaffected. At contact exit the pressure difference between the asperity point load and the EHD pressure spike of the line load will be smaller than the difference between the point load and the Hertzian line pressure on the entry side. Following, the local tensile surface stress at the asperity will be smaller at contact exit then it is at contact entry and the fatigue risk will therefore be relatively smaller at contact exit. Hence two important effects of lubrication on RCF; lubrication decreases the risk of competing damage modes and introduces asymmetry into the rolling contact pressure.

7. Conclusions The surface stress in front of a contact entering asperity was sufficient to initiate a RCF crack. After initiation the asperity point load mechanism could predict the spalling crack path. The crack grew at KI,max in the direction with the largest mode I crack tip load. Five crack path criteria were evaluated with identical results. It was concluded that the maximum mode I crack path was well approximated by the trajectory of the largest principal stress in the uncracked material when KI was the maximum during the load cycle. The simulated path agreed with the spalling profile both in the entry details and in the overall shape, which verifies the point load mechanism for RCF crack growth. Throughout the growth simulations, KI,max occured when the asperity entered the contact. With increasing length, the KII values of the mixed-mode cycle increased, which conformed to the irregular spalling profile and the wing cracks in the spall bottom. Qualitative arguments explained that a prospect kink down into the material was closed by negative T-stress and large negative mode I. A kink to the surface was not restricted by the same negative mode I load while the mode II load developed later in the cycle. Hence, the crack should eventually kink to form the spall. The RCF crack growth could be modelled with LEFM although the crack length was mostly small. The load conditions and mechanical material properties suggested that although the crack was short it was still larger than the estimated plastic zone at the crack tip and micro-structural dimensions of the pearlite steel. The simplified load cases with only normal asperity load, combined normal and tangential asperity load and the respective point force descriptions did not suffice to predict the spalling crack path but added understanding. The two point load

2866

D. Hannes, B. Alfredsson / Engineering Fracture Mechanics 78 (2011) 2848–2869

cases, with and without friction, could predict the two groups of surface distress crack paths. However, distributed pressures for both asperity and cylinder loads were needed for the correct path prediction. The gear contact could be simplified with a cylinder contact against a plane with an asperity. The FE method was required to separate and determine the evolvement of cylindrical and asperity pressures. When the contact pressures were known throughout the load sequence, the crack tip loads and the crack growth path in the spalling symmetry plane could be simulated numerically with elementary solutions. The elementary solutions gave high accuracy in the solutions while at the same time the computation time remained reasonable. The predictions for initiation and growth of RCF cracks followed accepted fatigue theory. The fatigue damage process was at both ends explained by the asperity point load mechanism and linear elastic fracture mechanics. Acknowledgement The authors greatfully acknowledge financial support from The Swedish Research Council. Appendix A. Closed-form stress fields The stress field in the symmetry plane was the superposition of elementary solutions for a half-plane with the cartesian coordinate system in Fig. 3. Observe that the cylindrical load was positioned at xd and a translation of the coordinate system was required. In the xz half-plane, a biaxial residual surface stress rR would be superimposed to rx. A.1. Cylindrical contact A normal concentrated line load, Pl (Flamant problem) at x = xd generates the following stress field [15]:

2P l

rx ¼  rz ¼ 

ðx  xd Þ2 z

p ððx  xd Þ2 þ z2 Þ2 2P l

sxz ¼ 

z3

p ððx  xd Þ2 þ z2 Þ2 2Pl

ðx  xd Þz2

p ððx  xd Þ2 þ z2 Þ2

;

ðA:1aÞ

;

ðA:1bÞ

:

ðA:1cÞ

McEwen (1949) [15] expresses the stress field below the Hertzian pressure in terms of m and n,

m2 ¼

qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi  f 2 þ 4ðx  xd Þ2 z2 þ f

ðA:2aÞ

qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi  f 2 þ 4ðx  xd Þ2 z2  f ;

ðA:2bÞ

1 2

and

n2 ¼

1 2

where f is

f ¼ a2l  ðx  xd Þ2 þ z2 :

ðA:2cÞ

The signs of m and n are the same as the signs of z and x  xd respectively. The stresses are

    p0l z2 þ n2  2z ; m 1þ 2 al m þ n2   p z2 þ n2 rz ¼  0l m 1  2 2 ; al m þn  2  p m  z2 sxz ¼  0l n 2 2 : al m þn

rx ¼ 

ðA:3aÞ ðA:3bÞ ðA:3cÞ

A.2. Asperity contact The stress field in the xz-plane from an inclined concentrated point force is given by the normal component, Pp (Boussinesqs problem) [15]:

D. Hannes, B. Alfredsson / Engineering Fracture Mechanics 78 (2011) 2848–2869

    Pp ð1  2mÞx2 z x2 z  3 ; 1  2p q4 q q5 3P z3 rz ¼  p 5 ; 2pq 3P xz2 sxz ¼  p 5 ; 2p q

rx ¼

2867

ðA:4aÞ ðA:4bÞ ðA:4cÞ

and the tangential component, Qp (Cerruti problem) [15]:

rx ¼

" # Qp x 3x3 3x x3 2x3 ð1  2mÞ 3  5  þ þ ; 2p q q qðq þ zÞ2 q3 ðq þ zÞ2 q2 ðq þ zÞ3

3Q p xz2 ; 2pq5 3Q p x2 z ¼ : 2pq5

ðA:5aÞ

rz ¼ 

ðA:5bÞ

sxz

ðA:5cÞ

Both stress fields are expressed in terms of

q¼

pffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi x2 þ z 2 :

ðA:6Þ

Explicit stress functions below a sliding spherical Hertz contact are given by Hamilton [16] and are expressed in terms of A, S, M, N, U, F and H defined as

A ¼ z2  a2p ; qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi S ¼ A2 þ 4z2 a2p ; rffiffiffiffiffiffiffiffiffiffiffiffi SþA ; M¼ 2 rffiffiffiffiffiffiffiffiffiffiffiffi SA ; N¼ 2 a  U ¼ arctan p ; M F ¼ M 2  N2 þ Mz  Nap ; H ¼ 2MN þ Map þ Nz:

ðA:7aÞ ðA:7bÞ ðA:7cÞ ðA:7dÞ ðA:7eÞ ðA:7fÞ ðA:7gÞ

For the normal Hertzian pressure distributions the stresses in the xz-plane are given by

rx ¼



 3Pp Mzap 1 2 ð1  2mÞ 3 ; ðNS þ 2AN þ a ð1 þ m Þz U  N  ð1  m ÞNz Þ  m Mza  p p x2 3 2pa3p S

  3Pp Mzap N þ ; 3 2pap S 

 3Pp xN xzH :  ¼ z S 2pa3p F 2 þ H2

ðA:8aÞ

rz ¼

ðA:8bÞ

sxz

ðA:8cÞ

These formulae are valid for x – 0. When x = 0,

" # a  a3p 3Pp p ; ð1 þ m Þðz arctan Þ þ  a p 2pa3p z 2ða2p þ z2 Þ " # a3p 3Pp ; ¼ 2pa3p a2p þ z2

rx ¼

ðA:9aÞ

rz

ðA:9bÞ

sxz ¼ 0:

ðA:9cÞ

The tangential Hertzian traction gives

rx

"

 2za3p 3Q p m Map Sm  2Am þ z2 x2 z2  m 2 ¼ x 1 þ U  ð1  2 m Þ þ  þ þ 1  x 2pa3p 4 3x3 x3 2 S 4 ( )# 2 2 2 2 2 x ap Nz ð1  2mÞðS þ 2AÞ z þ 3ap  ð7 þ mÞx þ þ ; þ 3 x 12 4 S

ðA:10aÞ

2868

D. Hannes, B. Alfredsson / Engineering Fracture Mechanics 78 (2011) 2848–2869

rz ¼

sxz

" ( )# x2 þ z2 þ a2p 3Q p Nz ; 1  2pa3p 2x3 S

" ( )#

3Q p 3zU 2 1 N z2  3a2p  x2  3ðS þ 2AÞ þ 2 : ¼ þ Mzap 2  x S x 2pa3p 2 4

ðA:10bÞ

ðA:10cÞ

Again these stress functions are only valid for x – 0. When x = 0,

rx ¼ 0; rz ¼ 0;

ðA:11bÞ

sxz ¼

ðA:11cÞ

" # a  3Q p 3 z2 ap p z arctan a þ  : p 2 2pa3p z 2ðz2 þ a2p Þ

ðA:11aÞ

Appendix B. Influence of total crack path curvature A through-thickness circular-arc edge crack in an infinite half-space was subjected to a remote uniform tension stress r0, see Fig. B.14a. The crack length was a and the crack front was assumed straight. The arc angle a gave a dimensionless measure of the total crack path curvature. For a through-thickness straight edge crack, a = 0. The normal stress,

rN ¼ r0 cos2 ðanÞ;

ðB:1Þ

and the tangential stress along the curved crack faces in an uncracked solid,

rT ¼ r0 cosðanÞ sinðanÞ;

ðB:2Þ

were used in Eq. (5) to compute the SIFs for a given a. The expressions for the SIFs were assuming a through-thickness straight edge crack, but were applied to curved edge cracks through the use of rN and rT. Dimensionless shape factors fI pffiffiffiffiffiffi and fII were obtained by normalizing the SIFs with r0 pa. Fig. B.14b presents the shape factors and thus the SIFs as function of a. For the applied load, increased total crack path curvature resulted in decreased KI and increased KII, which was consistent. Comparison with the shape factors proposed by Noda and Oda [35], showed fairly good agreement, especially for low values of a, see Fig. B.14b. For larger total crack path curvatures, KI computed by Eq. (5) gave an overestimation compared to KI by Noda and Oda. The results by Noda and Oda, fNO, were used to compute the absolute relative difference jf  fNOj/fNO, which was observed to increase with a. For a = 75, the difference for fI was approximately 35%. The opposite observation, but to a lower extent, was made for KII. For a = 75, the difference for fII was less than 9%. Appendix C. Influence of crack front curvature Newman and Raju [36] propose an empirical KI for a semi-elliptical surface crack, see Fig. C.15a. For a surface crack in an infinite half-space, subjected to a remote uniform tension stress r0 in the x-direction, KI at z = a is

Fig. B.14. Problem description and results illustrating the influence of total crack path curvature on the SIFs for a given applied load.

D. Hannes, B. Alfredsson / Engineering Fracture Mechanics 78 (2011) 2848–2869

2869

Fig. C.15. Problem description and results illustrating the influence of crack front curvature on KI.

pffiffiffiffiffiffi K I ¼ r0 paf ðuÞ;

1:13  0:09u with f ðuÞ ¼ pffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi : 1 þ 1:464u1:65

ðC:1Þ

The shape factor f was expressed in terms of u = a/c. Observe that u can be interpreted as a dimensionless measure of the crack front curvature at z = a. For a semi-circular crack u = 1, whereas for a straight crack u = 0. Note that the crack was straight in the xz-plane. Fig. C.15b shows how f and thus KI depend on u. Assuming that crack front curvature did not influence KI at z = a, resulted in an overestimation of KI when u was non-zero. The dashed line in Fig. C.15b represents the absolute relative difference with respect to a straight crack front, which reached 41% for a semi-circular surface crack. References [1] [2] [3] [4] [5] [6] [7] [8] [9] [10] [11] [12] [13] [14] [15] [16] [17] [18] [19] [20] [21] [22] [23] [24] [25] [26] [27] [28] [29] [30] [31] [32] [33] [34] [35] [36]

Dahlberg J, Alfredsson B. Int J Fatigue 2007;29:909–21. Tallian TE. Failure atlas for Hertz contact machine elements. New York: ASME Press; 1992. Way S. J Appl Mech 1935;2:A49–58. Ding Y, Rieger NF. Wear 2003;254:1307–17. Bastias PC, Hahn GT, Rubin CA, Gupta V, Leng X. Wear 1994;171:169–78. Bower AF. J Tribol 1988;110:704–11. Keer LM, Bryant MD. J Lubr Technol 1983;105:198–205. Bold P, Brown MW, Allen RJ. Fatigue Fract Engng Mater Struct 1992;15:965–77. Olsson M. In: Beynon J, Brown M, Smith R, Lindley T, Tomkins B, editors. Engineering against fatigue, proceedings of an international conference. The Netherlands, Sheffield: A.A. Balkema; 1999. Dahlberg J, Alfredsson B. Int J Fatigue 2008;30:1606–22. Alfredsson B, Olsson M. Engng Fract Mech 2000;65:89–106. Alfredsson B, Olsson M. Fatigue Fract Engng Mater Struct 1999;22:225–37. Hamrock BJ, Schmid SR, Jacobson BO. Fundamentals of fluid film lubrication. 2nd ed. New York: Marcel Dekker Corp.; 2004. Schneiders R, Schindler R, Weiler F. In: Proceedings of the 5th international meshing roundtable. p. 205–15. Johnson KL. Contact mechanics. Cambridge: Cambridge University Press; 1985. Hamilton GM. In: Proceedings of the institution of mechanical engineers, vol. 197C. Tada H, Paris PC, Irwin GR. The stress analysis of cracks handbook. Saint Louis, Missouri: Paris Production, Inc.; 1985. Qian J, Fatemi A. Engng Fract Mech 1996;55:969–90. Mróz K, Mróz Z. Engng Fract Mech 2010;77:1781–807. Erdogan F, Sih GC. J Basic Engng 1963;85:519–27. Sih GC. Mechanics of fracture 1 – Methods of analysis and solutions of crack problems. Noordhoff International; 1973. Palaniswamy K, Knauss WG. Int J Fracture 1972;8:114–7. Nuismer RJ. Int J Fracture 1975;11:245–50. Richard HA, Fulland M, Sander M. Fatigue Fract Engng Mater Struct 2005;28:3–12. Flodin A, Andersson S. Wear 1997;207:16–23. Suresh S, Ritchie RO. Int Metals Rev 1984;29:445–76. Riemelmoser FO, Pippan R. Int J Fracture 2002;118:251–70. Suresh S. Fatigue of materials. 2nd ed. Cambridge University Press, Massachusetts Institute of Technology; 1998. Rice JR. In: Fatigue crack propagation – ASTM STP 415. p. 247–311. Fett T. Engng Fract Mech 1997;57:365–73. Yu MT, Topper TH, Au P. In: Fatigue 84, proceedings of the 2nd conference on fatigue and fracture thresholds, vol. 1, Birmingham, UK. p. 179–90. Silva FS. Int J Fatigue 2004;26:241–52. Romeiro FFJ, Domingos CA, de Freitas MJM. In: McClung RC, Newman Jr JC, editors. Advances in fatigue crack closure measurement and analysis: second, vol. ASTM STP 1343. p. 321–36. Pommier S, Bompard P. Fatigue Fract Engng Mater Struct 2000;23:129–39. Noda N-A, Oda K. Int J Fracture 1993;64:239–49. Newman Jr JC, Raju IS. Engng Fract Mech 1981;15:185–92.