Numerical modeling of sub-surface initiated spalling in rolling contacts

Tribology International 59 (2013) 210–221

Contents lists available at SciVerse ScienceDirect

Tribology International journal homepage: www.elsevier.com/locate/triboint

Numerical modeling of sub-surface initiated spalling in rolling contacts Nick Weinzapfel n, Farshid Sadeghi Purdue University, School of Mechanical Engineering, West Lafayette, IN 47907, USA

a r t i c l e i n f o

a b s t r a c t

Article history: Received 5 August 2011 Received in revised form 30 December 2011 Accepted 14 March 2012 Available online 27 March 2012

A 3D finite element model was developed to investigate the influence of microstructure topology on the stochastic nature of rolling contact fatigue. Grains of the material microstructure are modeled with random Voronoi tessellations. Continuum damage mechanics and mesh partitioning are implemented to capture the initiation and propagation phases of fatigue damage that lead to spalling. Simulated fatigue spalling is shown to progress similarly to experimental observations of rolling contact fatigue. The fatigue lives obtained with the model exhibit scatter on par with empirical measures and are fit well by 2 and 3-parameter Weibull distributions. & 2012 Elsevier Ltd. All rights reserved.

Keywords: Rolling contact fatigue Finite element Damage mechanics

1. Introduction Material fatigue is the leading cause of failure in rolling element bearings (REBs) that are properly maintained and operate under conditions of elastohydrodynamic lubrication. As bearing elements roll between the raceway surfaces of the inner and outer rings, the ring material is fatigued by repeated exposure to non-conformal contact stresses. Hence, the phenomenon of material fatigue in REBs is commonly called rolling contact fatigue (RCF). It manifests through various mechanisms, but surface and subsurface initiated spalling are dominant [1,2]. Both mechanisms are characterized by fatigue cracks that propagate through the critically stressed volume to liberate material from the raceway surfaces [3] as shown in Fig. 1 [4]. However, with low frictional forces, clean lubrication, and good surface finishes, the operating conditions are favorable for subsurface initiated spalling. RCF is distinguished from classical fatigue in several ways that make it a challenging topic of research. The non-conformal contact between a rolling element and raceway produces a complex, multi-axial state of stress that diminishes as a function of distance from the contact site. Due to the rolling action, the directions of the principal stresses are not fixed for a given material point during the load cycle and the loading is nonproportional, i.e., the stress components do not rise and fall in

Abbreviations: AD, Anderson–Darling statistic; CVD, Carbon vacuum degassed; DEM, Discrete element model; FEM, Finite element model; FIB, Focused ion beam; JIC, Jump-in-cycles; LST, Linear strain tetrahedral; MP, Mesh partitioning; OSL, Observed significance level; RCF, Rolling contact fatigue; REB, Rolling element bearing n Corresponding author. E-mail addresses: [email protected] (N. Weinzapfel), [email protected] (F. Sadeghi). 0301-679X/$ - see front matter & 2012 Elsevier Ltd. All rights reserved. http://dx.doi.org/10.1016/j.triboint.2012.03.006

succession together. The localized nature of the contact precludes bulk failures and enhances the effects that heterogeneous microscale features such as grain size and distribution, inclusions, carbides, etc., have on the RCF process. Random variation of these inhomogeneities between otherwise identical bearings leads to scatter in their fatigue lives. Numerous empirical and research models have been developed to predict the fatigue lives of REBs [5]. Early efforts were focused on developing formulae to compute the ‘‘safe loads’’ and service lives of REBs from extensive full-scale fatigue tests [6,7]. Lundberg and Palmgren later formulated a more unified theory [8,9] that has since served as the basis for other models and ANSI/ ABMA and ISO standard load-life equations [10]. In general, fatigue life predictions of empirical models have their foundation in the stress solution of elastic contact and incorporate scatter through an assumed Weibull probability distribution function. On the other hand, research models provide deterministic life predictions by addressing the underlying physical mechanisms responsible for RCF. Typically, research models account for either the crack initiation or propagation phase but rarely both. An emerging class of numerical models of RCF provides an alternative to the conventional approaches that captures the stochastic nature of the phenomenon without relying on Weibull regression parameters obtained from full-scale bearing fatigue test data [11–14]. This new type of model utilizes randomly generated microstructure models and applies damage mechanics to simulate the fatigue induced material degradation. They capture both the initiation and propagation phases of RCF and have been used to obtain fatigue life estimates and spalling profiles that compare well with experimental results. These models have been formulated as discrete element (DEM) [11] and finite element (FEM) [12–14] models. The DEMs consider intergranular fatigue failure and treat the material grains as rigid

N. Weinzapfel, F. Sadeghi / Tribology International 59 (2013) 210–221

Nomenclature C D E F I N Rn T V Y a b dg e f k l m n p[max] tx, tz x xc

stiffness tensor [N/m2] damage tensor or scalar [dimensionless] modulus of elasticity [N/m2] failure probability [dimensionless] identity tensor [dimensionless] number of stress cycles elapsed [dimensionless] triaxiality function [dimensionless] traction vector acting on grain boundary [N/m2] element volume [m3] damage rate of strain energy release [N m] Hertzian contact half-width in the rolling direction [m] fatigue strength exponent [dimensionless] grain diameter [m] shape parameter of Weibull distribution [dimensionless] ratio of minimum strain energy release increment to total [dimensionless] sample size [dimensionless] Hertzian contact length transverse to the rolling direction [m] material parameter in damage rate law [dimensionless] grain boundary normal vector, number of elements [dimensionless] [maximum] pressure in the contact region [N/m2] shear, normal surface tractions [N/m2] surface coordinate where tractions are calculated [m] surface coordinate defining center of pressure distribution [m]

particles interconnected through a spring network [11]. The spring properties are calibrated such that the homogenized macro scale behavior simulates a continuum, and fatigue damage is propagated through the model by weakening and breaking the springs. Using the FEM approach, the material grains are deformable and damage may propagate along intergranular or transgranular pathways. A common simplification used in models for RCF is a 2D representation of the material. Though the contact loading may warrant such a simplification, the 3D internal geometry of the

z

211

depth below contact surface [m]

DN increment in elapsed stress cycles [dimensionless] Ds stress range [N/m2] Dtcritical critical shear stress reversal acting on a grain boundDjD

jD a b

e n s t tf0 m

ary [N/m2] increment in strain energy released by damage growth [N m] total strain energy released by damage growth [N m] location parameter of Weibull distribution [dimensionless] scale parameter of Weibull distribution [dimensionless] strain tensor [dimensionless] Poisson’s ratio [dimensionless] stress tensor or scalar [N/m2] shear stress [N/m2] shear fatigue strength coefficient [dimensionless] coefficient of friction [dimensionless]

Subscripts H hydrostatic c contact e, i, j, k, l, m, n indices eq equivalent (von Mises) r fatigue resistance s surface D portion of DN attributable to required damage increment j portion of DN attributable to required strain energy release increment

material should be considered for its role in the formation and propagation of fatigue damage. Spall surface topologies,.e.g., Fig. 1, and reconstructions of sub-surface fatigue cracks from focused ion beam (FIB) images [15] are evidence of this process. This work presents a numerical model of RCF developed in a 3D finite element framework that explicitly considers the microstructure topology for its stochastic influence on the fatigue of REBs. Grains of the material microstructure are modeled using randomly generated Voronoi tessellations, and damage mechanics is incorporated to account for the gradual material degradation due to the cyclic contact fatigue loading. The model uses a mesh partitioning procedure to simulate the initiation, propagation, and coalescence of fatigue damage in the form of intergranular cracks. Simulation times are significantly expedited by an integration algorithm that considers the strain energy released through damage formation and growth. The spalling patterns generated by the model for rolling line contact loading originate under the raceway and propagate to the surface with three-dimensional characteristics that are unobtainable with plane strain/stress models. Examination of sequential crack patterns reveals the influence of the 3D microstructure topology on the failure progression. The topological effect on fatigue lives is discussed and the results are compared with empirical observations.

2. Simulation of a contact stress cycle in a roller bearing

Fig. 1. Typical fatigue spall in AISI 52100 (from [4]).

The rolling contact fatigue model presented in this investigation was developed to simulate subsurface initiated spalling in cylindrical roller bearings. The normal contact between a roller and a raceway can be modeled as an equivalent contact between

212

N. Weinzapfel, F. Sadeghi / Tribology International 59 (2013) 210–221

two cylinders, which produces a pressure distribution given by p ¼ pmax

rffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi  x 2 1 a

ð1Þ

that acts over a rectangular region of length, l, and half-width, a. To compute the subsurface stress history in the raceway of a cylindrical roller bearing, a finite element model of a 3D elastic half-space, depicted in Fig. 2, is constructed using the microstructure model generator developed in [16]. The dimensions are limited to (10  1  12)a in the rolling (X), transverse (Y), and vertical (Z) directions, respectively. All displacements on the lower boundary of the model are fixed, and periodicity constraints are enforced in the rolling and transverse directions in order to effectively simulate the response of a semi-infinite, elastic half-space. The topology of the material microstructure is modeled by randomly created Voronoi tessellations. A Voronoi tessellation is produced by placing a set of generating points in space and dividing the volume into cells with the characteristic that all points in a particular cell are closer to the generating point of that cell that any other member of the generating set [17]. This procedure yields set of convex polyhedral shapes that represents the material grains to a good degree of accuracy [18–20]. The grain size, dg, is controlled by specifying a minimum distance criterion between closest neighbor points of the generating set. Randomly placing the generating points ensures that the assembly of grains is different for each microstructure model. Grains are discretized into tetrahedra by placing a node at the centroid of each face to divide the faces into triangles forming the bases of tetrahedra that share a common apex at the grain centroid as shown in Fig. 3. Each tetrahedron is meshed with

Fig. 3. (a) Voronoi cell model of material grain, (b) tetrahedral discretization.

either linear or quadratic elements depending on its proximity to the applied loads as described in [16]. The microstructure model is subjected to a sequence of surface tractions to simulate the action of a bearing element rolling over the raceway. The pressure distribution of Eq. (1) is applied at a sequence of 21 discrete locations from  2a to 2a as depicted in Fig. 2 using the equation sffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi   ðxxc Þ 2 t z ðxÞ ¼ pmax 1 ð2Þ a where x is the coordinate where the traction is evaluated, and xc is the location of the center of the pressure. Varying xc enables the traction distribution to translate over the surface of the halfspace. Shear surface traction in the rolling direction is defined by a constant ratio of the normal traction: t x ðxÞ ¼ ms 9t z ðxÞ9

ð3Þ

where m is the coefficient of friction. The model’s response to these tractions is evaluated using ABAQUS/Standard [21].

3. Fatigue damage model Repeated exposure to rolling contact stresses results in a progressive deterioration of the material through the nucleation of growth of micro-voids and cracks. Damage mechanics provides a convenient framework for modeling the influence of these microscopic failure mechanisms on the homogenized response of a material. The material degradation can be coupled to stress– strain constitutive equations through a thermodynamic state variable for damage, D. In general, the constitutive relationships take the following form:

sij ¼ C ijkl ðIklmn Dklmn Þemn

ð4Þ

where sij, Cijkl, Iklmn, Dklmn, emn are the stress, stiffness, identity, damage, and strain tensors, respectively. For the case of isotropic damage, the damage tensor reduces to a scalar, simplifying Eq. (4) to

sij ¼ C ijkl ð1DÞekl

ð5Þ

where 0 rD r 1

Fig. 2. Finite element model of contact stresses applied to an elastic half-space.

ð6Þ

with a value of D ¼0 corresponding to a pristine material and D ¼1 for the fully damaged state. The evolution of damage at a material point must be characterized by a rate law appropriate for the phenomenon of interest. Bearing steels that fail by high-cycle, rolling contact fatigue exhibit quasi-brittle behavior for which the following

N. Weinzapfel, F. Sadeghi / Tribology International 59 (2013) 210–221

damage rate law is appropriate [22]  m dD Ds ¼ dN sr ð1DÞ

213

ð7Þ

where Ds is the stress range that promotes the failure, and sr and m are material parameters. The quantity, sr, is called the resistance stress, because it characterizes a material’s ability to resist accumulation of fatigue damage [23]. The stress range Ds of interest is chosen based on the material’s response to the contact loading. In the absence of strong shear tractions acting at the contact site, mode I crack growth is suppressed by the compressive hydrostatic component of the stress field. Therefore, it is believed that shear stress reversals are ultimately responsible for RCF [3,8,9,24–26]. The present work assumes that the grain boundaries constitute weak planes in the material, and the action of reversing shear stress on these planes advances the fatigue failure. After simulating a rolling contact stress cycle, the elemental Cartesian stresses, s, acting at the grain boundaries are expressed with respect to their local grain boundary coordinate systems, s0 . The traction vector, T0 , acting on the grain boundary is then determined to be 8 0 9 s = > < 11 > 0 ð8Þ T 0 ¼ s0 Un0 ¼ s12 > ; : s0 > 13

0

where n is the normal vector of the grain boundary surface. The two shear components are resolved into a single shear vector that traces a path during the rolling contact cycle as exemplified in Fig. 4. Using the longest chord method [27], the critical shear stress reversal, Dtcritical, is determined and applied in Eq. (7) to compute the damage evolution rates of the elements as   dD Dtcritical m ¼ ð9Þ dN tr ð1DÞ The critical shear stress reversal, Dtcritical, depends on the location of the material point, as in the usual solution of elastic contact, as well as the orientation of the grain boundary plane

Fig. 5. Torsion S–N data for through-hardened (58-62 HRC), JIS SUJ2 from [28] with power law fit.

passing through the point. The collection of grain boundaries varies among the different microstructure models so that fatigue damage evolves differently in each domain. Hence, the dispersion of fatigue lives manifests naturally from the explicit representation of the microstructure topology with Voronoi tessellations. The two material parameters, tr and m, in Eq. (9) are calibrated using test data from torsion fatigue which is presumed to accumulate damage by shearing mechanisms similar to RCF. In this investigation, torsion S–N data for through-hardened bearing steel JIS SUJ2 (AISI 52100 equivalent), shown in Fig. 5 (from [28]), is used to extract the material parameters as follows. A powerlaw fit of the data has the form of Basquin’s rule

Dt=2 ¼ t0f Nb ¼ 2:3879N 0:090092

ð10Þ

where t is the fatigue strength coefficient, and b is the fatigue strength exponent. Integrating (9) from D ¼0 to D ¼1 yields 1  tr m tr -Dt ¼ N1=m ð11Þ N¼ m þ 1 Dt ðm þ 1Þ1=m 0 f

Comparing Eqs. (10) and (11) yields 1 ¼ 11:1 b   1 b tr ¼ 2t0f 1 ¼ 5979 MPa b

m¼

ð12Þ

4. Micro crack initiation and propagation

Fig. 4. Example tip trace of resolved shear stress vector on grain boundary plane illustrating measurement of Dtcritical.

When elements reach a threshold level of damage, intergranular fatigue cracks are inserted in the microstructure model using a newly developed mesh partitioning (MP) procedure. The MP procedure creates physical discontinuities, i.e., crack surfaces, (where grains debond) in the model by duplicating the nodes shared by neighboring elements in the mesh while leaving the shape of the mesh intact. Before describing the method, the concepts of physical entities and the relationships between them are defined. The physical entities considered by the MP procedure are the geometric features of the grains, i.e., faces, edges, and vertices. Faces are composed of the element surfaces on the periphery of a grain, and the unique sets of edges and vertices collected from the faces constitute the edges and vertices of the grain. Each face, edge, and vertex is associated with its parent grain and retains a unique identity in the microstructure assembly. Initially, the congruent entities of neighboring grains deform and displace equally. As fatigue damage evolves and cracks are formed, the relationships between the congruent entities change

214

N. Weinzapfel, F. Sadeghi / Tribology International 59 (2013) 210–221

and allow them to behave independently. Three abstract concepts, joints, chains, and clusters, are defined to characterize and track the statuses of the relationships between faces, edges, and vertices of neighboring grains, respectively. A joint defines the relationship between a pair of faces of two neighboring grains. The faces deform identically under load while the joint is intact. When the joint’s status is set to broken, a crack is introduced between two grains, and the faces are allowed to deform uniquely. Two types of joints that must be considered: (1) internal joints defined between congruent faces of adjacent grain pairs; and (2) external joints assumed to exist between faces that lie on the exterior of the microstructure model. The latter type is broken by default but makes possible a unified and efficient treatment of chains. Chains characterize the relationship(s) between two or more congruent edges of adjacent grains. Initially, all the edges in a chain deform as a single edge; however, they may deform independently after the joints linking their parent grains are broken. Chains are categorized as internal, if the edges traverse through the interior of the model, and external, if the edges lie on the exterior surface(s). A cluster characterizes the relationship(s) among two or more vertices of adjacent grains that coincide at a single point in the non-partitioned mesh. The displacement of such a group of vertices is initially defined by the degrees of freedom of a single node in the finite element mesh. As the joints are broken between the vertices’ parent grains, the vertices involved in the cluster may become disassociated from one another, thus requiring multiple nodes to describe their displacements. The MP procedure is implemented in a three-step approach that sequentially examines the statuses of joints, chains, and clusters to determine when conditions warrant duplication of the nodes. The three steps of MP are demonstrated by cutting in half a block specimen composed of eight cubic grains shown in Fig. 6. Fig. 6(a) and (b) illustrates the arrangement of the cubes and its associated finite element mesh of linear strain tetrahedra (LST). In order to cut the block, the four joints in Fig. 6(c) are broken. The first step examines the status of each internal joint to determine if a discontinuity should be incorporated into the mesh to model a crack. If the joint is declared broken, the procedure illustrated in Fig. 7 is used to introduce a crack between two grains and decouple the deformations of the faces. First, the nodes that are unique to the joint are identified. The nodes on the perimeter of the joint are not unique to it, because they may be involved in other joints. Copies are made of the unique nodes and assigned to grain 2, leaving the original nodes associated with the face on grain 1. Application of step 1 to the four broken joints shown in Fig. 6(c) produces crack surfaces in the finite element model as illustrated in a cut view of the block in Fig. 6(d). External joints need not be considered because they are merely abstractions and do not couple the deformations of any grain faces. Step 2 examines the statuses of the joints in each chain to determine which, if any, edges or groups of edges should be separated from the others in the chain. A chain must contain at least two broken joints, internal or external, to be eligible for consideration in step 2. Fig. 6(e) illustrates that four internal chains (blue) and eight external chains (green) are affected by the broken joints in the cubic grain example. For each internal chain, there are four distinct grain edges that must be separated into two pairs of edges as a result of cutting the block in half. Examining the list of edges and joints corresponding to an internal chain (Table 1) determines the number times to duplicate the original node. The joints and edges are ordered in the table as they would be encountered. For examples: joint 1 is between edges 1 and 2; joint 4 is between edges 4 and 1. Each joint’s status is written above it: ‘‘I’’¼ intact, ‘‘B’’ ¼broken. The first edge

Fig. 6. Mesh partitioning illustrated for a block specimen: (a) arrangement of cubic grains, (b) meshed grains, (c) broken joints, (d) cut view of specimen after step 1 showing internal crack surfaces, (e) chains affected by broken joints, (f) specimen after step 2, (g) clusters affected by broken joints, (h) specimen after step 3, (i) nodes duplicated.

encountered in the chain retains the original node by default. Joint 1 is broken, therefore edge 2 and any following edges, prior to encountering another broken joint, are assigned copy 1 of the node. Any edge appearing after the last broken joint defaults to using the original node. In the block example with LST elements,

N. Weinzapfel, F. Sadeghi / Tribology International 59 (2013) 210–221

215

Table 2 Joints and related vertices for the central cluster of Fig. 6(g). Joint

1

2

3

4

5

6

7

8

9

10

11

12

Status Vertex 1 Vertex 2

I 1 5

I 2 6

I 3 7

I 4 8

B 1 2

B 5 6

B 4 3

B 8 7

I 1 4

I 2 3

I 5 8

I 6 7

Table 3 Step 3 of MP — identifying vertex groups within central cluster of Fig. 6(g). Vertex Group Group Group Group Group Group Group Group

Fig. 7. Step 1 of the mesh partitioning procedure.

Table 1 Geometric and relational information for an internal chain. Status Joint list Edge list Node copy

B

I

1 1 0

B

2 2 1

I

3 3 1

4 4 0

there is only one node between the endpoints of each edge; therefore, only one new node per edge pair required. Elements with higher order interpolation functions would require additional new nodes. Fig. 6(f) depicts the specimen in tension after step 2. The grain faces and edges have been decoupled; only the connections between grain vertices remain to be addressed. Fig. 6(g) illustrates the sites of clusters affected by the four broken joints in Fig. 6(c). Step 3 of the MP procedure determines which nodes should be duplicated to decouple the displacements of vertices on neighboring grains. This is accomplished by assuming that all of the vertices are uncoupled and examining the status of the joints between the parent grains to determine which vertices can be connected. The central cluster in Fig. 6(g), with its eight vertices (one per grain) and twelve joints, is used to demonstrate this process. Table 2 describes the set of joints, their statuses, and associated pairs of vertices on the parent grains bordering each joint for this cluster. Table 3 demonstrates the logic used to identify groups of vertices that are constrained together. Each vertex is assumed to be in its own group to start, then the vertices associated with each intact joint are grouped together.

1 – – – – – – – –

start joints 1–4 joints 5–8 joint 9 joint 10 joint 11 joint 12 end

2

3

4

5

6

7

8

0 1 2 3 4 5 6 7 0 1 2 3 0 1 2 3 These joints are broken — no change in groupings 0 1 2 0 0 1 2 0 0 1 1 0 0 1 1 0 Vertices 5 and 8 in same group already — no change Vertices 6 and 7 in same group already — no change 0 1 1 0 0 1 1 0

Joint 1 is intact and connects vertices 1 and 5 (see Table 2) which are in groups 0 and 4, respectively. The lower of the two group numbers is then assigned to both of these vertices, such that vertex 5 belongs to group 0. Additionally, any other vertex whose current group assignment matches the larger group number will be assigned the smaller group number. Joints 1–4 each regroup a single vertex. Joints 5–8 are broken and therefore do not change the vertex groupings. The intact status of joint 9 regroups vertex 4, which is currently in group 3, into group 0. Notice, that vertex 8 is also currently in group 3, so it too joins group 0. Joint 10 has similar effect, while joints 11 and 12 enforce redundant constraints. Finally, group 0 contains vertices 1, 4, 5, and 8 that are associated with the original node, while group 1 contains vertices 2, 3, 6, and 7 which are associated with a copied node. Application of step 3 to each of the clusters in Fig. 6(g) decouples the displacements of the grain vertices on opposite sides of the cut and frees the halves of the block to move independently as depicted in Fig. 6(h). The nodes duplicated in each step of MP are identified in Fig. 6(i). After intergranular cracking occurs, grains are prevented from interpenetrating and their sliding behavior is characterized by a simple Coulomb friction model. This is configured in ABAQUS using surface-to-surface CONTACT PAIRS with a desired coefficient of friction, mc, and a HARD PRESSURE-OVERCLOSURE relationship.

5. Numerical implementation The rolling contact fatigue simulation involves simultaneously solving the damage evolution equation (9) and the constitutive relationships (5) for the microstructure model. Evaluation of Eq. (9) requires knowledge of the critical shear stress reversal, Dtcritical, which would be determined ideally for every stress cycle. Such an approach would be impractical, because rolling element bearings can endure millions or billions of stress cycles; therefore, a compromise must be made. A standard ‘‘jump-incycles’’ (JIC) algorithm [29] has been used in 2D simulations of RCF [11–14], wherein the stress-strain response of the model is assumed to be constant over a finite number of cycles, DNi, in a block i. A constant increment in damage, DD, is typically prescribed for the element experiencing the largest damage rate per cycle in each block, thereby governing the number of cycles in the

216

N. Weinzapfel, F. Sadeghi / Tribology International 59 (2013) 210–221

block i, DNi. Examination of Eq. (9) reveals that the rate of damage per cycle approaches infinity as D approaches 1. This behavior can cause DNi to be less than 1, making the standard JIC algorithm slow to progress. In this section a modified JIC algorithm is presented to avoid the situation just described. The modification is motivated by two simple ideas: (1) strain energy is released by the formation and growth of damage as fatigue progresses, (2) the strain energy released in each block of cycles should advance the total by some minimum fraction of the previous total. The modified JIC procedure is performed as follows: (i) The number of elapsed cycles and released strain energy is set to 0 at the beginning of the RCF simulation. Also, the material is assumed to be in its pristine state, so that damage in all elements is 0. Ni ¼ 0 jD ¼ 0 Die ¼ 0,

e ¼ 1,. . .,nelements

ð13Þ

(ii) Critical shear stress reversals Dtie,critical are computed for each element for the current block of cycles. (iii) Damage evolution rates for each element can then be determined. #m  i " dD Dtie ¼ ð14Þ dN e tr ð1Die Þ (iv) Additionally, the maximum strain energy release rates for the elements are computed by Y ie ¼

ð ieq,e Þ2 Rin,e V e 2Eð1Die Þ2

s

where Ve is the element volume and  2 2 sH Rin,e ¼ ð1 þ nÞ þ3ð12vÞ 3 seq

ð15Þ

ð16Þ

with sH and seq measured at the instant during the contact stress cycle where seq is maximum. (v) The element experiencing the maximum damage evolution rate is chosen as the most critical for the current block of cycles. "  #  i i dD dD ¼ max ð17Þ dN critical dN e

supplemented by a number of cycles, DN if . 9 8 0, DjiD =jD Zf > > > > = < ! i f jD DN f ¼ i i > 1 DN D , DjD =jD o f > > > ; : Dji D

DN i ¼ DNiD þ DNif

ð20Þ

(ix) The damage states and strain energy release increment are updated as required by (viii).  i dD Dieþ 1 ¼ Dieþ 1 þ DNif dN e  i dD DjiD ¼ DjiD þ SY ie DNif ð21Þ dN e The additional number of cycles is applied in small subincrements so that when elements reach the critically damaged state, D¼1, they are excluded from subsequent computations. (x) If any element has reached the critically damaged state during this cycle block, then its corresponding joint is broken and a fatigue crack is introduced in the microstructure using the MP procedure. That is if Dieþ 1 Z 1-mcrack ðbreak the associated jointÞ

ð22Þ

and the states of damage for elements associated with the broken joint are set to 0 and cannot change thereafter. (xi) Finally, the total released strain energy and the number of elapsed cycles are updated.

jD ¼ jD þ DjiD N ¼ N þ DNi

ð23Þ

A value for f must be chosen to define the minimum increment of strain energy released relative to the total released in previous cycle blocks. The modified JIC algorithm is identical to the standard method when f is set to 0, but causes large errors in the solution if f is chosen to be arbitrarily large. Increments of strain energy released during each block of cycles using the standard JIC algorithm are shown in Fig. 8. The percentage change can often exceed 2% of the accumulated total, but remains close to 1% using a 5-point moving average. Therefore, a choice was made to set f¼2% for the RCF simulation in order to complement, but

(vi) The number of cycles for the current block is computed as 

DNiD ¼ DD=

dD dN

i ð18Þ critical

(vii) The state of damage in each element is updated and the strain energy release during this block of cycles is calculated. i dD DN iD dN e  i dD DjiD ¼ SY ie DNiD dN e

Dieþ 1 ¼ Die þ



ð19Þ

(viii) The increment in strain energy release is checked against a desired minimum increment relative to the total strain energy released thus far. If it is not at least as large as the threshold, then the number of cycles in the block is

Fig. 8. Strain energy release increment as a percentage of the accumulated strain energy released for the standard jump-in-cycles algorithm.

N. Weinzapfel, F. Sadeghi / Tribology International 59 (2013) 210–221

not consistently exceed the strain energy released by the standard JIC algorithm. Application of the modified JIC algorithm substantially reduces the number of computational cycles required to complete the analysis without causing significant deviations from the results obtained with the standard JIC algorithm. Fig. 9 illustrates that the trends in the number of elapsed stress cycles, critically damaged elements, and total released strain energy remain nearly identical despite the significant reduction in the number of computational cycles. The sharp split in Fig. 9(d) indicates the point in the simulation where the modified JIC algorithm took effect. Prior to the split, the two approaches behaved identically because the strain energy increment predicted by the standard procedure surpassed the 2% minimum requirement. The crack patterns produced using both approaches are very similar as shown in Fig. 10. Table 4 summarizes the computational benefits for the same microstructure domain. The reduction in processor (CPU) time depends non-linearly on the number of computational cycles

217

saved. This is because convergence of the finite element analysis occurs at a slower rate in later computational cycles as the number of cracks increases. The computational expense is reduced by a factor of nearly 5 for the domain used in the comparison. Investigating the role of microstructure topology on fatigue life scatter requires simulations with several domains. Accomplishing this feat becomes much more practical by using the modified JIC algorithm.

6. Results and discussions This section describes the results obtained from the 3D FEM for sub-surface initiated spalling in cylindrical roller bearings. The localized nature of the non-conformal contacts experienced in REBs precludes bulk failure of the bearing rings and enhances the influence of heterogeneous microscale features on the fatigue process. Typical dimensions of the contact region between the bearing element and raceway can be on the order a hundred to a few hundred microns. At the lower end of the spectrum, these dimensions are comparable to characteristic lengths of subsurface topology features in bearing steels where grain sizes are on the order of 10 mm. The same grain size is used to randomly construct the 33 microstructure models studied in this investigation. Table 5 summarizes the simulation parameters. The surface tractions of Eqs. (2) and (3) are prescribed to the model with a maximum Hertzian pressure of 2.0 GPa and friction coefficient of 0.05 that is representative of lubricated rolling contacts. The material properties are assumed to be isotropic and homogeneous so that the only differences between the 33 domains are found in the distribution of the material grains and their geometric features, or topological disorder. 6.1. Spall pattern

Fig. 9. Comparison of modified and standard jump-in-cycles algorithms.

Fig. 10. Comparison of crack patterns generated using the standard and modified jump-in-cycles algorithms.

Table 4 Computational expense for standard and modified JIC algorithms. Approach

Standard JIC Modified JIC

Computational cycles

316 121

Clock time (h)

297 63.24

CPUs

8 8

CPU time (h)

(days)

2376 506

99.0 21.1

Fig. 11 illustrates the progression of failed joints in a sample microstructure model as a function of increasing rolling contact cycles. The images show the model cut at three different planes along the transverse rolling direction and rotated 201 toward the viewer so that the upper surface is partially visible. Each column depicts the evolution of cracks for a particular cross-section, and each row shows the variation in the crack patterns through the transverse rolling direction at a particular stress cycle count. The first crack appears below the surface near the depth where the maximum orthogonal shear stress occurs (z/a ¼0.5) consistent with the observations of low-friction, rolling contact [24,25]. Notice that it is visible in only one cross-section due to the finite dimensions of the material grains, whereas a crack in a 2D model (plane strain condition) implies instantaneous failure across all cross-sections in the direction transverse to rolling. As the number of stress cycles increases to 244 million, more isolated cracks appear and short range linking occurs. The accumulation of damage continues in this fashion producing long horizontal Table 5 Rolling contact fatigue simulation parameters. Material grain diameter (dg) Elastic modulus (E) Poisson’s ratio (n) Hertzian line contact half-width (a¼ 5 dg) Maximum Hertzian contact pressure (pmax) Surface coefficient of friction (ms) Subsurface coefficient of friction in cracks (mc) Fatigue resistance stress, tr Damage rate law exponent, m Damage increment (DD) Minimum increase in released strain energy, f

10mm 200 GPa 0.3 50 mm 2.0 GPa 0.05 0.4 5979 MPa 11.1 0.2 2%

218

N. Weinzapfel, F. Sadeghi / Tribology International 59 (2013) 210–221

Fig. 11. Evolution of sub-surface crack patterns at various transverse planes due to rolling contact fatigue.

cracks below the surface that propagate toward the surface in some instances. Although the cracks may surface in some zcross-sections, the 3D topology inhibits the cracks from instantly breaching the surface across the entire width of the microstructure. After 715 million cycles, the sub-surface damage has saturated and the surface cracks have grown nearly wide enough to produce a spall formation. Crack propagation drawings for GCr15 bearing steel tested in rolling contact fatigue experiments [24] depict behavior that qualitatively agrees well with the model. 6.2. Effect of microstructure topology on fatigue lives To study the influence of subsurface topological disorder on the fatigue lives, the 33 microstructure domains were analyzed with homogeneous material properties, i.e., uniform E, n, tr and m. Identical simulation parameters, given in Table 5, were used in the fatigue simulation of each microstructure model. Although the stress-strain responses of all the homogeneous, pristine microstructure models are the same when expressed in the XYZ reference frame, there is significant variation in the stresses acting on the weak planes or grain boundaries. This variation appears in the critical shear stress reversal Dtcritical experienced by each element, causing non-uniform evolution of the constitutive

relationships via Eqs. (5) and (9). Therefore, the progression from a pristine material to fatigue failure is a unique process for each microstructure model. For the purposes of this study, the number of cycles elapsed prior to the first crack’s appearance constitutes the initiation life for a domain. The additional number of cycles experienced up to the formation of a surface crack is called the propagation life, and its sum with the initiation life is the final life. The depth below the surface where the first fatigue crack appears in each domain is shown in Fig. 12. The influence of the topological disorder can be seen in the variation of this quantity. Two points of interest are (1) that the average crack initiation depth, z/a¼ 0.529, is very close to the depth of the maximum orthogonal shear (0.5) stress as calculated with Hertz theory and (2) that the range of locations varies within the bounds of observations made by Chen et al. [25] (0.329rz/ar0.7). Fig. 13 illustrates the lives of the 33 domains decomposed into their initiation and propagation phases. It can be seen that the lives are typically dominated by the propagation phase, and the initiation phase constitutes an average 31% of the total life. The initiation lives show little variation, and the scatter in the final lives are primarily governed by the propagation phases. A number of cumulative distribution functions (CDFs) could be used to characterize the scatter of fatigue lives of REBs; however, the form most commonly applied is the family of Weibull CDFs

N. Weinzapfel, F. Sadeghi / Tribology International 59 (2013) 210–221

219

Fig. 12. Depth of crack initiation in the rolling contact fatigue model. Fig. 14. Weibull probability plot of initiation and final lives with 2 and 3-parameter Weibull fits.

Fig. 15. Relative magnitude and directions of critical shear stress reversals plotted on the unit sphere according to the orientation of the corresponding grain boundary: (a) single domain, (b) composite from several domains. Fig. 13. Initiation and propagation phases of the microstructure models.

[30]. The general form of the Weibull CDF is given by, 

 Na e , N 4a FðN; a, b,eÞ ¼ 1exp 

b

ð24Þ

where F is the probability that failure occurs by a number of cycles, N, while a r0, b o0, eo0 are the location, scale, and shape parameters, respectively, that define the distribution function. The parameter a describes the finite life that all members of a population can expect to survive with 100% reliability. The minimum life parameter is typically excluded (a ¼0) in analyses of fatigue lives of REBs, reducing Eq. (24) to the popular 2-parameter Weibull distribution. For the 2-parameter Weibull distribution, e is called the Weibull slope and measures the slope of the line Fig. 14    1 ¼ elnðNÞelnðbÞ ln ln ð25Þ 1F that relates life cycles to failure probability. Therefore, e provides some measure of the scatter in the fatigue life data; the smaller the slope, the broader the scatter. Fig. 15 displays the initiation and final lives on a Weibull probability plot with 2 and 3-parameter fits. The Weibull slopes for the initiation and final lives were determined to be 22.42 and 3.88, respectively. Fatigue data for REBs manufactured from carbon vacuum-degassed (CVD) AISI

52100 steel follow a range of Weibull slopes between 0.51 and 5.7 [31] which bounds the model’s predictions for final lives. Endurance testing of nine angular contact ball bearings made from CVD AISI 52100 steels with a maximum contact pressure of 3.1 GPa attained lives between 1e8 and 1e9 stress cycles before forming spalls 1–2 mm in size [4]. Possible reasons for the model attaining similar lives with a 2.0 GPa maximum contact pressure include the following: (1) plasticity is not considered in the model and therefore, residual stresses have been neglected (2) the analysis was stopped once the first crack reached the surface, rather than waiting until a similar sized spall was formed. A comparison of the results generated by the current model and previous computational models of RCF is presented in Table 6. The value of L10 ¼243 million cycles is higher than calculated by earlier computational models primarily due to two reasons. The first is that the damage rate law parameters, tr and m, of the current model are based on torsion S–N data for modern bearing steel [28]. The lives of these specimens are nearly an order of magnitude larger than the 1951 data [32] used to calibrate the modeling parameters in [11–13]. Note that the L10 life of the current model also exceeds the results of models [14] calibrated with the modern torsion S–N data. This is attributable to the fact that the grain boundaries in 2D models are generally exposed to larger stress reversals than in the current 3D model. Therefore, the damage propagation occurs more slowly in the 3D model. Fig. 15 helps elucidate the reason for this effect by

220

N. Weinzapfel, F. Sadeghi / Tribology International 59 (2013) 210–221

Table 6 Comparison of results obtained with computational models of RCF (CZ ¼cohesive zone). Reference

Model type

k domains

sr (MPa)

m

einit

efinal

L10 (  1e6 cycles)

Average % initiation

[11] [12] [13] [14] [14] Current model

2D 2D 2D 2D 2D 3D

40 40 35 40 40 33

6113 6113 6113 5979 5979 5979

10.1 10.1 10.1 11.1 11.1 11.1

11.75 5.11 7.15 2.93 2.94 22.42

1.85 4.08 12.54 2.74 3.14 3.88

23.47 49.42 16.12 39.1 37.7 243

19.8 (10–20)a 85.9 30 33 31

a

explicit DEM implicit FEM explicit FEM explicit FEM with linear CZ explicit FEM with bilinear CZ implicit FEM with MP

Value was estimated from source though not reported within.

illustrating the magnitude of the critical shear stress reversals as they relate to the orientations of the grain boundaries where they occur in the pristine microstructure. The shear stress reversals are largest on the grain boundary planes parallel to the X and Z planes and least on those parallel to the Y plane. In the 2D models, grain boundaries cannot have orientations of the latter type and are restricted to have normal directions that lie on the Y plane as indicated by the arc in Fig. 15(b). Therefore, a grain boundary in a 2D model is more likely to experience a larger shear stress reversal than one in the current 3D model, thus the propagation of damage to the surface is faster for a 2D model. The initiation lives are less scattered than the final lives as is generally the case for [11,12,14], and the fraction of the life spent in the initiation phase is similar to all but [13]. The likely explanation for this is that the states of fully damaged finite elements were allowed to remain at D ¼1 in [13]. Kotzalas noted several published studies in which RCF failures correlate well with 2-parameter Weibull distributions over midrange probabilities but deviate toward a minimum fatigue life expectancy at high levels of reliability [33]. The lives of the microstructure models also demonstrate this behavior suggesting the use of a 3-parameter Weibull distribution to fit the data. The parameters for the final lives were obtained using the method of maximum likelihood estimate (MLE) [34]: a ¼1.04e8 cycles, b ¼3.24e8 cycles, and e¼2.77. The goodness-of-fit for the 3-parameter distribution can be assessed with the Anderson–Darling statistic (AD) AD ¼ k þ

k X

  1 ð12iÞ ln½FðNi Þ þ ln 1FðN k þ 1i Þ ¼ 0:380 ki¼1

After adjusting for the sample size (k¼33)   0:2 AD0 ¼ AD 1þ pffiffiffi ¼ 0:393 k

ð26Þ

computational expense is necessary to solve the damage-coupled, elasticity equations hundreds of times during a simulation. This motivated the development of a modified ‘‘jump-in-cycles’’ algorithm to expedite the solution which reduced the computational effort by a factor of nearly 5 for the case demonstrated. Results generated with 33 different microstructure domains agree well both qualitatively and quantitatively with experimentally obtained data. Failure was shown to originate at various depths below the surface due to the random microstructure topology. After a relatively short initiation phase (31% of the total life), isolated cracks propagated and linked to form long cracks beneath the surface in line with experimental observations. Eventually the cracks breached the surface, but development of realistic spall pattern was delayed by the finite 3D geometric features of the microstructure. The scatter in RCF lives of the domains as measured by the Weibull slope was found to be within the limits recorded for carbon vacuum degassed AISI 52100. The final life data appeared to shown a tendency to deviate from a 2-parameter Weibull regression toward a finite life at high levels of reliability as similarly shown in empirical studies. Using a 3-parameter Weibull distribution proved to fit the data well.

Acknowledgments The authors would like to express their thanks to the sponsors of the Mechanical Engineering Tribology Laboratory for the generous support and insightful discussions that made this work possible. References

ð27Þ

the observed significance level (OSL) is determined to be

  OSL ¼ 1= 1þ exp 0:1 þ1:24lnðAD0 Þ þ 4:48AD0 ¼ 0:376 ð28Þ This is much greater than the critical value of 0.05 needed to accept the hypothesis that the data comes from a Weibull distribution with parameters listed previously, thereby confirming the adequacy of the 3-parameter Weibull fit.

7. Summary and conclusions This paper presented a finite element model for investigating the influence of 3D microstructure topology on the stochastic phenomenon of rolling contact fatigue. Randomly generated Voronoi tessellations were used to represent the topological disorder present in bearing steels at the micro scale level. The model incorporated damage mechanics to simulate the progressive failure of JIS SUJ2 (AISI 52100 equivalent) subjected to cyclic contact fatigue loading. The initiation and propagation of fatigue cracks (grain debonding) was considered with a mesh partitioning procedure. Significant

[1] Littmann WE, Widner RL. Propagation of contact fatigue from surface and sub-surface origins. ASME Journal of Basic Engineering 1966;88:624–36. [2] Littmann WE. The mechanism of contact fatigue. NASA Special Report 1969;SP 237. [3] Tallian TE. Failure Atlas for Hertz Contact Machine Elements. New York: ASME Press; 1999 pp. 179-231. [4] Rosado L, Forster NH, Thompson KL, Cooke JW. Rolling contact fatigue life and spall propagation of AISI M50, M50NiL, and AISI 52100, Part I: experimental results. STLE Tribology Transactions 2010;53:29–41. [5] Sadeghi F, Jalalahmadi B, Slack T, Raje N, Arakere NK. A review of rolling contact fatigue. Journal of Tribology 2009;131(4):041403. [6] Goodman J. Roller and ball bearings. Proceedings of the Institute of Civil Engineering 1912;189:82–166. [7] Palmgren A. Die Lebensdauer von Kugellagern. Zeitschrift des Vereines Deutscher Ingenieure 1924;68(14):339–41 (see also The service life of ball bearings. NASA Technical Translation, TT 1-13460; 1971). [8] Lungberg G, Palmgren A. Dynamic capacity of rolling bearings. Acta Polytechnica, Mechanical Engineering Series 1947;1(3):1–50. [9] Lundberg G, Palmgren A. Dynamic capacity of roller bearings. Acta Polytechnica, Mechanical Engineering Series 1952;2(4):1–32. [10] Barnsby R, Duchowski J, Harris T, Ioannides S, Loesche T, Nixon H, Webster M. Life Ratings for Modern Rolling Bearings — A Design Guide for the Application of International Standard ISO 281/2. New York: ASME International; 2003 TRIB-Vol. 14. [11] Raje N, Sadeghi F, Rateick RG. A statistical damage mechanics model for subsurface initiated spalling in rolling contacts. Journal of Tribology 2008;130(4):042201.

N. Weinzapfel, F. Sadeghi / Tribology International 59 (2013) 210–221

[12] Jalalahmadi B, Sadeghi FA, Voronoi FE. Fatigue damage model for life scatter in rolling contacts. ASME Journal of Tribology 2010;132(2):021404. [13] Slack T, Sadeghi F. Explicit finite element modeling of subsurface initiated spalling in rolling contacts. Tribology International 2010;43:1693–702. [14] Slack T, Sadeghi F. Cohesive zone modeling of intergranular fatigue damage in rolling contacts. Tribology International 2011;44:797–804. [15] Grabulov A, Ziese U, Zanbergen HW. TEM/SEM investigation of microstructural changes within the white etching area under rolling contact fatigue and 3-D crack reconstruction by focused ion beam. Scripta Materialia 2007;57: 635–8. [16] Weinzapfel N, Sadeghi F, Bakolas V. An approach for modeling material grain structure in investigations of Hertzian subsurface stresses and rolling contact fatigue. Journal of Tribology 2010;132(4):041404. [17] Okabe A, Boots B, Sugihara K, Chiu SN. Spatial Tessellations: Concepts and Applications of Voronoi Diagrams.2nd edn. West Sussex: John Wiley & Sons; 2000. [18] Ito O, Fuller ER. Computer modeling of anisotropic grain microstructure in two dimensions. Acta Metallurgica et Materialia 1993;41(1):191–8. [19] Zavattieri PD, Espinosa HD. Grain level analysis of crack initiation and propagation in brittle materials. Acta Materialia 2001;49:4291–311. ¨ [20] Mucklich F, Osher J, Schneider G. Die Charakterisierung homogener poly¨ ¨ edrischer Gefuge mit hilfe des Raumlichen Poisson-Voronoi-Mosaiks und der Vergleich zur DIN 50 601 (The characterization of homogeneous polyhedral microstructures applying the spatial Poisson-Voronoi tessellation compared ¨ Metallkunde 1997;88(1):27–32. to the standard DIN 50 601). Zeitschrift fur [21] Dassault Syste mes. 2009, Abaqus 6.9 Documentation, Dassault Syste mes Simulia Corp, Providence, RI, USA. http://abaqusdocs.ecn.purdue.edu:2080/ v6.9/. [22] Xiao YC, Li S, Gao Z. A continuum damage mechanis model for high cycle fatigue. International Journal of Fatigue 1998;20(7):503–8.

221

[23] Bolotin VV, Belousov IL. Early fatigue crack growth as the damage accumulation process. Probabilistic Engineering Mechanics 2001;16(4):279–87. [24] Chen L, Chen Q, Shao E. Study on initiation and propagation angles of subsurface cracks in GCr15 bearing steel under rolling contact. Wear 1989;133: 205–13. [25] Chen L, Shao E, Zhao D, Guo J, Fan Z. Measurement of the critical size of inclusions initiating contact fatigue cracks and its application in bearing steel. Wear 1991;147:285–94. [26] Zaretsky EV, Parker RJ, Anderson WJ. A study of residual stress induced during rolling. ASME Journal of Lubrication Technology 1969;91:314–9. [27] Lemaitre J, Chaboche JL. Mechanics of Solid Materials. Cambridge: Cambridge University Press; 1990. [28] Shimizu S, Tsuchiya K, Tosha K. Probabilistic stress-life (P–S–N) study on bearing steel using alternating torsion life test. Tribology Transactions 2009;52:807–16. [29] Lemaitre J. A Course on Damage Mechanics. Berlin: Springer-Verlag; 1992. [30] Weibull W. A statistical theory of the strength of materials, Ingeniors Vetenkaps Akademien. Proceedings of the Royal Swedish Academy of Engineering 1939;151: 1–45. [31] Harris TA, Barnsby RM. Life ratings for ball and roller bearings. Proceedings of the ImechE, Part J: Journal of Engineering Tribology 2001;215(6):577–95. [32] Styri H. Fatigue strength of ball bearings races and heat-treated 52100 steel specimens. Proceedings of the American Society for Testing and Materilas 1951;51:682–700. [33] Kotzalas MN. Statistical distribution of tapered roller bearing fatigue lives at high levels of reliability. Journal of Tribology 2005;127:865–70. [34] Lockhart RA, Stephens MA. Estimation and tests of fit for the three-parameter Weibull distribution. Journal of the Royal Statistical Society, Series B (Methodological) 1994;56(3):491–500.