A discrete damage mechanics model for high cycle fatigue in polycrystalline materials subject to rolling contact

International Journal of Fatigue 31 (2009) 346–360

Contents lists available at ScienceDirect

International Journal of Fatigue journal homepage: www.elsevier.com/locate/ijfatigue

A discrete damage mechanics model for high cycle fatigue in polycrystalline materials subject to rolling contact Nihar Raje, Trevor Slack, Farshid Sadeghi * School of Mechanical Engineering, Purdue University, 585 Purdue Mall ME, West Lafayette, Indiana 47906, USA

a r t i c l e

i n f o

Article history: Received 11 June 2008 Received in revised form 4 August 2008 Accepted 8 August 2008 Available online 24 August 2008 Keywords: Damage mechanics Discrete material model Fatigue life scatter Polycrystalline materials

a b s t r a c t Fatigue behavior of polycrystalline materials is significantly influenced by their microstructural topology. The microstructural heterogeneity is one of the primary reasons for dispersion in high cycle fatigue lives of such materials. In this work, a damage mechanics based fatigue model that incorporates gradual material degradation under cyclic loading is presented in conjunction with a discrete material representation that takes the material microstructural topology into account. Microstructures are generated stochastically through the process of Voronoi tessellation. Micro-crack initiation, coalescence and propagation stages are modeled using damaged zones in a unified framework. The model is applied to study high cycle fatigue in rolling contacts. The effect of material topological disorder and inhomogeneity on fatigue life dispersion is studied. Fatigue damage is found to originate sub-surface and propagate towards the surface. Sub-surface damage patterns from the model are consistent with experimental observations. Propagation life is found to constitute a significant fraction of total life. Lives are found to follow a 3parameter Weibull distribution. The relative proportion of lives spent in the initiation and propagation stages are in good quantitative agreement with experiments. Ó 2008 Elsevier Ltd. All rights reserved.

1. Introduction Commercially used bearing steels exhibit polycrystalline microstructures that are made up of a large number of individual crystals or grains that are separated by grain boundaries. The stochastic and heterogeneous nature of such microstructures due to different grain structures creates an inherent topological disorder in the material that can lead to significant scatter in macro material characteristics, e.g. fracture and fatigue responses. The modeling of such topological disorder can be conveniently done in the framework of discrete material models [1,2]. Discrete material models have been extensively used to study fracture processes in brittle materials [3–7]. Such models do not require assumptions on the initial size and distribution of microcracks as in continuum mechanical models based on fracture mechanics. Rather, the nucleation, coalescence and growth stages of cracks are modeled in a unified framework as an outcome of simulations. The two main types of discrete models available in literature are the spring network or lattice models [8–11,7] and the cohesive zone models [3,4,12]. In lattice models the material is represented by a lattice of network sites or nodes that are connected by simple mechanical elements such as springs or beams. Failure is incorpo* Corresponding author. E-mail addresses: [email protected] (N. Raje), [email protected] (T. Slack), [email protected] (F. Sadeghi). 0142-1123/$ - see front matter Ó 2008 Elsevier Ltd. All rights reserved. doi:10.1016/j.ijfatigue.2008.08.006

rated through the breaking of these elements under a specified stress or energy criteria. Material disorder is modeled through the assignment of individual material properties to each lattice element according to a statistical distribution [13,9,14] or by using random lattice geometry [11,15]. Cohesive zone models on the other hand are based on a phenomenological framework in which the material fracture characteristics are described by a cohesive surface traction-displacement relationship in the process zone. Failure due to fracture is achieved by a progressive softening of material in the process zone under a prescribed cohesive law. The use of discrete models to study fatigue crack initiation and growth under cyclic loading has received considerable attention in the last decade [16–20]. However, most of the efforts have been directed towards deterministic modeling of the fatigue phenomena using cohesive zones in the framework of the finite element method. Such models albeit providing a detailed description of the failure process are computationally expensive. Recently, Rinaldi et al. [2] proposed a lattice model to simulate polycrystalline microstructures and studied the resultant scatter in fatigue response at the macro level. Their model offered significant computational advantages over full scale finite element models based on a discretization of microstructure at the sub-grain level. However, the model involved the use of one dimensional linear normal elastic springs as lattice elements representing grain boundaries. Thus, the shear component of force that could lead to grain boundary sliding was ignored. Further, the addition of tangential springs leads to more accurate force distributions in the lattice. Also, their

N. Raje et al. / International Journal of Fatigue 31 (2009) 346–360

analysis involved the use of a simplified Basquin type stress-life relationship for the springs and linear damage accumulation according to the Miner’s rule. Experimental results for damage evolution in metals under cyclic loading suggest a highly non-linear behavior with respect to number of load cycles [21]. Therefore, the study of fatigue requires the use of more sophisticated damage evolution models. This paper presents a stochastic fatigue model that combines the ideas of discrete material representation and continuum damage mechanics to study the initiation and propagation of fatigue damage in polycrystalline materials. Note that the model is limited to intergranular crack initiation and propagation. The model incorporates non-linear fatigue damage accumulation and is conceptually along the lines of the spring lattice models. However, interfaces connecting lattice sites are modeled using continuous normal and tangential spring sets instead of one-dimensional normal springs. This gives normal, tangential and torsional stiffnesses to these interfaces. Material topological disorder is incorporated in the model through the use of Voronoi tessellations to simulate the polycrystalline material microstructure. In addition, material property variation is modeled through the assignment of individual properties to each interface in accordance with a statistical distribution. Two sources of fatigue life scatter are considered: (1) the topological randomness due to geometric variability in the material microstructure and (2) material property randomness due to a statistical distribution of material properties spatially throughout the material. 2. Material model The model is based on a discrete representation of the polycrystalline material [22] where the macro material domain is assumed to be formed by an assemblage of discrete, interacting, microelements that represent the grain geometries, see Fig. 1. The idea is to simulate only the geometric disorder due to the grain structure. Thus, the elements themselves are rigid but are connected along their contacting sides through fictitious fibers to form compliant joints (interfaces) as shown in Fig. 2. These fibers remain in their undeformed positions when there is no relative motion between elements, but undergo a compression or elongation otherwise. The joint between the contacting elements is assumed to be formed by infinite such fibers along its length. The fibers are modeled as visco-elastic members having certain stiffness and damping characteristics. In addition, for generality, the stiffness of these fibers is assumed to be different in the normal (n) and tangential (t) directions, where the two directions are as shown in Fig. 2. Note that the n direction is perpendicular to the contacting side of the element under consideration and points outwards. The

347

normal and tangential stiffnesses of the fibers are denoted as Kn and Kt respectively. The corresponding damping ratios are denoted as nn and nt respectively. Note that the normal and tangential damping coefficients are related to the stiffnesses by;

pffiffiffiffiffiffiffiffiffiffi C n ¼ nn 2 K n m;

pffiffiffiffiffiffiffiffiffi C t ¼ nt 2 K t m

ð1Þ

where, m is the average mass per unit length of the discrete elements. The normal and tangential spring stiffnesses, Kn and Kt, are functions of the overall elastic properties of the solid, namely, the Young’s modulus E and the Poisson’s ratio m. Potapov et al. [23] proposed the following relations expressing Kn and Kt, in terms of E and m for Voronoi elements;

Kn ¼

E ; Lp ð1  mÞ

Kt ¼

Eð1  3mÞ Lp ð1  m2 Þ

ð2Þ

Here, Lp is the average length of inter-element joints in the domain. Fig. 3a shows a typical interaction occurring between a pair of discrete elements. The two initially adjacent elements i and j, undergo a relative displacement with respect to each other due to application of external forces on the boundary. Thus, initially coincident points A (on element i) and B (on element j) are displaced with respect to each other under relative motion of the elements. Note that the relative displacements are shown highly exaggerated for illustrative purposes and are in practice small compared to the size of the micro-elements. Consider the fiber AB which is stretched by dn in the n direction and dt in the t direction at a given time step. The total displacement between points A and B is;

^ þ dt ^t d~ s ¼ dn n

ð3Þ

^ and ^t are the unit vectors along the n and t directions, where, n respectively. Since the elements are assumed to be rigid, the displacements dn and dt vary linearly along the length of the contacting side. The corresponding spring force generated in the fiber is given by;

^ þ K t dt ^t d~ F s ¼ K n dn n

ð4Þ

Note that the normal and tangential spring forces act in the same direction as the corresponding displacements dn and dt. This force acts on element i at point A as shown in Fig. 3b. An equal and opposite force acts on element j at point B. Note that, in general, the spring force will not act along the fiber because the spring stiffnesses in the normal and tangential direction are different. The damping force in the fiber is proportional to the time rate of change of the fiber displacement and is given by;

dn dt ^ þ C t ^t n d~ F d ¼ Cn ds ds

ð5Þ

where, dn and ddts are the time rates of change of the fiber displaceds ments in the normal and tangential direction, respectively. The damping force opposes relative motion between the elements and hence acts in the direction opposite to the displacement rate vectors. Fig. 3c shows the damping force acting on element i at point A. An equal and opposite force acts on element j at point B. Note the amount of damping is small and the force in the fibers is mainly due to the spring effect. The resultant force on element i at point A, located at a distance x from the midpoint M of the contacting side (see Fig. 4a) is;

d~ F ¼ d~ F s þ d~ Fd

ð6Þ

The total force acting on the contacting side of element i is then found by integrating the normal and tangential forces in each of the connecting fibers. It is given by;

~ F¼ Fig. 1. Discrete material representation.

Z

   Z L    dn dt ^þ K n dn þ C n K t dt þ C t dx n dx ^t ds ds L L L

ð7Þ

348

N. Raje et al. / International Journal of Fatigue 31 (2009) 346–360

Fig. 2. (a) Inter-element contact in the discrete model; (b) fiber model.

where, 2L is the length of the contacting side. The resultant force is assumed to act through the midpoint M of the side. Note that this is accompanied by a moment about M as shown in Fig. 4b, given by;

  dn ~M ¼ xdx M K n dn þ C n ds L Z

L

ð8Þ

The force ~ F has normal and tangential components Fn and Ft, respectively. The normal and tangential stresses on the joint are approximated by dividing the corresponding force components by the contact length and are shown in Fig. 5. Thus,

rn ¼

Fn Ft ;s ¼ 2L 2L

ð9Þ

The effect of the moment term is neglected since it is small in magnitude [22]. The joint forces and moments are then transferred to an equivalent force and moment at the centroid G of the elements associated with the joint, see Fig. 6. The computation is done for each element by summing the contributions of all joints associated with the element. The equivalent force and moment acting at the element centroid are given by;

~ FG ¼

n X j¼1

~ Fj;

~G ¼ M

n h X

~M ð~ r MGj  ~ FjÞ þ M j

i

ð10Þ

j¼1

where, ~ F j is the total force acting at the midpoint of a joint associated with the element, n is the number of joints on the element and ~ rMGj is the vector from G to the midpoint M of the joint. The equivalent force and moment are then used in the dynamic equilibrium equations for the element as illustrated in Section 5.

3. Generation of material microstructure Solid polycrystalline metals are formed through grain growth during the process of crystallization from their molten state. If during this process all material grains grow simultaneously and at the same rate then the resulting grain structure is well described by Voronoi tessellation [24]. In fact, the geometry of polycrystalline microstructures can be simulated to a good degree of accuracy through the process of Voronoi tessellation [25,26]. The process consists of using a set of randomly placed points as nucleation points or seeds and constructing regions around them such that all points enclosed by the region are closer to the given nucleation point than any other nucleation point. The resultant regions are convex polygons known as Voronoi polygons [27]. Due to the random distribution of the nucleation points, non-unique material domains are obtained for every tessellation simulation. This gives rise to a topological randomness in the microstructure. In the present model, the micro-elements (grains) of the discrete material domain are generated through such a Voronoi tessellation process. The density of nucleation points is defined in accordance with the grain size of the polycrystalline material. The generated element shapes have variable number of sides and variable orientations. However, the most probable number of sides is six, which also corresponds to the number of sides for maximum thermodynamic stability of material grains. Care is taken during the tessellation process so that the element sizes are uniform, i.e., there is not much variation in the grain size across the domain. This is achieved by setting

349

N. Raje et al. / International Journal of Fatigue 31 (2009) 346–360

Fig. 4. (a) Resultant fiber force, (b) Resultant joint force and moment.

σn

Element i

τ

Element j Fig. 5. Stress components on an inter-element joint.

Fig. 3. (a) Fiber displacement, (b) Fiber spring force, (c) Fiber damping force.

upper and lower limits on the distances between the nucleation points. 4. Damage mechanics model Fatigue damage is manifested through a progressive degradation of the material due to microscopic mechanisms like formation and growth of micro-cracks and voids. Damage mechanics [28,29] provides a convenient framework to treat these microscopic failure mechanisms in an empirical fashion. It leads to a better understanding of rupture problems in structures by the definition of an internal damage variable D at a material point which represents the gradual deterioration of materials before the initiation of a crack. The initiation and early growth of damage is discontinuous and is strongly affected by the heterogeneous nature of polycrystalline materials [30]. The progressive damage accumulation is

Fig. 6. The resultant force and moment on the center of gravity of the element.

350

N. Raje et al. / International Journal of Fatigue 31 (2009) 346–360

effected by introducing the damage variable into the material constitutive stress-strain relationship. For the most general elastic case, the damage modified constitutive equation takes the form,

rij ¼ C ijkl ðIklop  Dklop Þeop

ð11Þ

where, rij is the stress tensor, Cijkl is the material stiffness tensor containing the elastic constants, Dklop is a fourth order damage tensor and eop is the strain tensor. For a one-dimensional case, Eq. (11) reduces to,

r ¼ ð1  DÞEe

ð12Þ

Thus, the presence of damage leads to a corresponding stiffness reduction for the material. The damage variable D takes values ranging from 0 (undamaged) to 1 (completely damaged). Thus,

06D61

ð13Þ

A value of D = 1 corresponds to a state of micro-crack initiation at the concerned material point. In high cycle fatigue, the evolution of damage at a material point can be related to the stress level at that point through a non-linear equation of the form,

dD ¼ f ðr; DÞ dN

ð14Þ

where, N is the number of stress cycles and r is a stress measure. Note that Eq. (14) ignores the effects of microplasticity produced by local stress concentrations that originate from variability in geometry and properties at the microscopic level. For a more accurate analysis, the local plastic behavior of the material should be taken into account. This is beyond the scope of the present work. A widely used form of Eq. (14) for one-dimensional fatigue derived from thermodynamics is [31–34],

dD ¼ dN



m Dr rr ð1  DÞ

ðK t Þ ¼ ðK t Þ0 ð1  DÞ

dD ¼ dN



Ds rr ð1  DÞ

m ð17Þ

Here, Ds is the shear stress range acting along the inter-element joint. Since the cracks are always subjected to normal compressive loading the crack closure parameter h is set to 0 in this case. 5. Numerical implementation The stress and displacement fields in the discontinuous material domain are obtained through the solution of the dynamic equations of motion (EOM) for each discrete element. Analysis is performed in the framework of the Discrete Element Method (DEM) [37–39,23] where the EOM are solved for each element independently instead of solving a global set of matrix equations. The solution is affected through an explicit time integration scheme. In addition to the local relative damping between elements, a mass proportional global damping is assumed. The global damping coefficient Cg is given by;

C ig ¼ a  mi

ð18Þ

where, i = 1, 2 and 3 correspond to the x, y and h directions, respectively and a is a constant. The element velocities and positions are found using an explicit central difference scheme as;

x_ itþ1 ¼ 2

! ! 1  a D2t i Dt F it _ þ x t12 Dt Dt mi 1þa 2 1þa 2

xitþ1 ¼ xit þ Dt  x_ itþ1 2

ð15Þ

where, Dr is the critical stress range, rr and m are temperature dependent material parameters that have to be empirically determined. Note, rr which is in general a function of the mean stress is also referred to as the resistance stress [34] because it controls the resistance to damage accumulation. The damage variable is implemented within the current modeling framework through the reduction of the inter-element fiber stiffness components as

ðK n Þ ¼ ðK n Þ0 ð1  hDÞ

shear [36]. Thus, the parameter rr is assumed to be independent of the mean shear stress. Accordingly, an equation of the following form is used for fatigue damage evolution under rolling contact loading:

ð16Þ

where (Kn)0 and (Kt)0 are the initial normal and tangential joint stiffnesses under pristine conditions as illustrated in Fig. 7. h is a parameter that accounts for crack closure. The crack closure parameter is introduced to capture the effects of the partial closing of micro-cracks that occur in compression. Thus, h = 1 in tension and 0 6 h < 1 in compression. A value of h = 0 corresponds to complete closure of microcracks. Note that damage is assumed to be isotropic and hence the damage variable is scalar, i.e., the same damage quantity D is used for the reduction of normal and tangential stiffnesses. The critical stress quantity Dr responsible for fatigue damage is chosen based on the type of loading experienced by the structure. In the case of rolling contacts, the sub-surface normal stresses on all inter-element joints in the material domain are compressive and hence inhibit crack opening [35]. Therefore, it is assumed that the normal compressive stress acting on the joints does not cause any fatigue damage and the only stress component that favors the formation and propagation of cracks is the range of shear stress acting along the joints. This is along the lines of the idea that cracks are only propagated by mode II loading (mode I absent). It has also been observed that the mean stress does not have any effect in

ð19Þ ð20Þ

For stability purposes, the time step Dt is kept less than the natural period of free oscillation of the spring-mass system formed by the elements and connecting fibers [39]. Thus,

1 1 6 Dt 2p

rffiffiffiffiffiffi Kn m

ð21Þ

Based on the new element positions and velocities, the displacement and rates of displacement for the inter-element fibers are calculated for the next time-step. Inter-element joint stresses are then computed according to Eq. (9). Fatigue analysis involves the simultaneous solution of the damage evolution Eq. (17), stiffness degradation Eq. (16) and constitutive equations for each inter-element joint to account for the coupling between damage and the material constitutive behavior. Thus, the stress-strain relationships need to be solved and the joint stiffnesses need to be updated after every stress cycle to incorporate the current state of damage in the joint. However, it is impractical to perform the tremendous number of increments needed for the large number of stress cycles arising in high cycle fatigue. Therefore, a method that involves the use of the ‘jump-in-cycles’ procedure [21] is used for achieving computational efficiency. The method assumes a piecewise periodic, constant amplitude loading, i.e., the stress histories for the inter-element joints are assumed to remain unaltered over a finite number of cycles DN constituting a block i. The increment in damage DD is assumed to be constant over a block of cycles. Thus, the damage evolution is assumed to be piecewise linear with respect to the number of cycles (see Fig. 8). It is to be noted that the shape of the damage evolution curve is not predetermined but it is an outcome of the numerical simulations due to the stress-damage coupling. The value of DD is a tradeoff between accuracy and computation time. For higher resolution, DD needs to be chosen as small as possible, however

351

N. Raje et al. / International Journal of Fatigue 31 (2009) 346–360

σn

0

D = 0, (K n) = (K n) N

0

0 < D < 1 , (Kn ) = (K n )(1 -h D) (K n)

Δσ n

ε

D = 1, (K n ) = 0 Crack closure effect

τ

0

D = 0, (Kt ) = (Kt ) N

0

0 < D < 1 , (K t ) = (K t )(1- D) (K t )

Δτ

γ

D = 1, (K t ) = 0

Fig. 7. Degradation of joint stiffnesses with damage accumulation, (a) normal direction and (b) tangential direction.

(i) Initial damage states for each joint are assigned. For a pristine material domain containing no initial flaws, the initial damage in each joint is assumed to be 0. Thus,

1

D0j ¼ 0; j ¼ 1::::::::njoints

ð22Þ

The number of cycles elapsed is initialized to 0. Thus,

N¼0

D

(ii) Stress histories Dr for each joint are computed for the current block i. (iii) The damage evolution rate in each joint is evaluated knowing the stress history for the present block and the current state of damage in the joint using,

Δσ i

Δσ 0

0

ΔN

0

Δσ 1 ΔN 1



ΔN i

ð23Þ i j

dD dN

N

Joint Life Fig. 8. Piecewise linear approximation for damage evolution.

it directly affects the computation time. Stress calculations are performed once for every block of cycles and not for each individual cycle N. This saves a considerable amount of computation time since the major part of the computational expense is spent in solving for the stress fields in the domain. The number of cycles in a block is variable and it is obtained using the following procedure:

"

i ¼ j

Drij

#m

rr ð1  hDij Þ

ð24Þ

(iv) The joint with the maximum damage evolution rate is chosen as the critical joint for the current block of cycles. Thus,



dD dN

i crit

    dD i    ¼ Max   dN j 

ð25Þ

(v) The number of cycles in the current block of cycles is computed as,

DD DN i ¼  i dD dN crit

ð26Þ

352

N. Raje et al. / International Journal of Fatigue 31 (2009) 346–360

(vi) The number of cycles elapsed is updated to,

N ¼ N þ DN

i

ð27Þ

(vii) The increment in damage for each joint during the current block of cycles is then given by,

DDij ¼



i dD DN i dN j

ð28Þ

(viii) The damage states for each joint at the start of the next block of cycles are updated to,

Diþ1 ¼ Dij þ DDij j

ð29Þ

(ix) The joint stiffnesses are modified at the start of the next block of cycles according to, i

ðK n Þij ¼ ðK n Þ0j ð1  hDj Þ

ð30Þ

ðK t Þij ¼ ðK t Þ0j ð1  Dij Þ

Steps (i) through (ix) are successively repeated for each block of cycles. A micro-crack is assumed to be initiated at a joint when the accumulated damage in the joint equals 1. For subsequent calculations, the damage variable in the joint is set to 1, i.e.

Dij ¼ 1; if Dij > 1

ð31Þ

Note that this corresponds to a complete loss of tangential stiffness for the joint and hence the inability to carry tangential stress henceforth. The joint also loses the ability to support normal tensile loads. However, the joint can still support compressive loads through crack closure. At this stage, the normal joint stiffness in compression is reduced to

ðK n Þij ¼ ðK n Þ0j ð1  hÞ

ð32Þ

and is left unaltered for subsequent computations. The failed joint constitutes a damaged zone. The crack propagation phase is modeled through link up of damaged zones that correspond to multiple micro-cracks coalescing together to form longer cracks. 6. Evaluation of damage parameters The fatigue damage evolution equation introduces two new material parameters rr and m that have to be experimentally determined. The method used here identifies these parameters from experimental stress-life data (S-N curve) available in literature. The general equation for an S-N curve under completely reversed one-dimensional loading takes the form of Basquin’s law given by,

 B Dr 2A ¼ AðNf Þ1=B ) Nf ¼ 2 Dr

ð33Þ

Here, Nf is the number of cycles to failure at the stress level Dr while A and B are material constants. An integration of Eq. (15) assuming the stress level to be unchanged and neglecting crack closure effects gives,

Z 0

Nf

dN ¼

Z 0

1



m

rr ð1  DÞ Dr

dD ) Nf ¼

h r im 1 r ðm þ 1Þ Dr

ð34Þ

Comparison of Eqs. (33) and (34) gives,

m ¼ B; rr ¼ 2AðB þ 1Þ1=B

ð35Þ

Analysis for rolling contact is performed using AISI-52100 through hardened bearing steel which has a polycrystalline microstructure

and is extensively used in rolling element bearings. In order to evaluate the damage parameters in this case, equivalence between torsional fatigue and rolling contact fatigue is assumed. It is hypothesized that since fatigue damage under rolling contact conditions is caused purely by the action of shear stresses the mechanism of failure is similar to that prevalent in torsional fatigue. Thus, by knowing the stress-life relation (S-N curve) for a material in torsional fatigue, the fatigue damage parameters can be extracted. From the data for completely reversed torsional fatigue for AISI-52100 steel [40], the following parameters are obtained:

m ¼ 10:1; rr ¼ 6113 MPa

ð36Þ

7. Model application This section presents the application of the proposed fatigue model to rolling contacts where the scale of the domain is comparable to the microstructural length scale and the topological disorder due to different grain structures can significantly affect the fatigue life. Rolling contacts arise when two non-conformal bodies roll with respect to each other producing alternating stressing over very small volumes. They are experienced in applications such as rolling element bearings, gears, cam-followers and rail-wheel contacts. The alternate stressing makes the contacting bodies susceptible to failure due to fatigue, a phenomenon commonly referred to as rolling contact fatigue [41,42]. Typical dimensions over which the bodies make contact in such situations range from 100– 500 lm. Furthermore, the stress fields under such loading are highly localized around the contact region. As a result, the fatigue failure risk is highly localized in the volume surrounding the contact region. The dimensions of the volume susceptible to fatigue are comparable to the microstructural dimensions, e.g., typical grain sizes for AISI-52100 bearing steel used in rolling bearing applications are of the order of 10 lm. This is the average microelement size used in analysis. Since the model is two-dimensional, it is applied to a rolling line contact (plane strain condition) which typically arises in roller bearings, gears and cam-follower systems. In non-conformal contacts, the contact width (2b) is small as compared to the radii of curvature of the contacting bodies. Under such conditions, the stresses in the contact region are neither critically dependent on the shape of the contacting bodies distant from the contact area, nor upon the precise manner in which they are supported [43]. The stresses can be calculated to a good approximation by considering each body as semi-infinite in extent and having a plane surface on which the appropriate contact pressure acts. In accordance with this, one of the contacting bodies is replaced by a semi-infinite domain as shown in Fig. 9 and its interaction with the other body is modeled by applying an appropriate pressure distribution on the contact surface of the domain. In theory, a semiinfinite domain is infinite in extent. However, for computational purposes, the domain created is restricted to 10b in the x (rolling) direction and 7b in the z (depth) direction. Here 2b is the Hertzian contact width. The limits are selected in consideration of the fact that in Hertzian contacts, stresses at points sufficiently remote from the contact region are practically unaffected by the contact loading. Fixed boundary conditions are imposed on the bottom surface of the semi-infinite domain in order to maintain equilibrium of the micro-elements. The side boundaries are left free. Simulation parameters used in this investigation are given in Table 1. Analyses are carried out for bearing steel AISI-52100 using 40 material domains each generated stochastically using the process of Voronoi tessellation. For each domain, the number of elements is kept constant. This number is determined based on the average distance between nucleation points and the overall domain size. A surface traction coefficient of l = 0.05 that is representative of typical lubricated contacts is chosen. All analyses are carried out using

N. Raje et al. / International Journal of Fatigue 31 (2009) 346–360

a Hertzian pressure distribution on the surface with a corresponding shear traction q related to the normal pressure p by,

Table 1 Simulation parameters for rolling contact study

qðxÞ ¼ l  pðxÞ

Hertzian contact half-width (b) Elastic modulus of material (E) Mean micro-element size Number of micro-elements in domain Number of joints in domain Damage increment (DD) Surface friction coefficient (l)

ð37Þ

Note that the pressure is applied dynamically using an exponential rise function. The loading history experienced by points in the subsurface region of a rolling contact can be simulated to a good degree of accuracy by moving the pressure distribution along the top surface of the semi-infinite domain as shown in Fig. 9.

353

100 lm 200 GPa 10 lm 4500 13885–13914 0.1, 0.2 0.05

7.1. Mechanism of failure Failure in rolling contacts is manifested through the formation of pits or spalls on the surface of contacting bodies. To study the process through which this occurs, analysis is carried out using constant elastic properties throughout the domain for a maximum

Hertzian pressure of pmax = 2 GPa. Fig. 10a shows the development of damaged zones (failed inter-element joints) with increasing number of rolling contact cycles in a sample material domain. It is seen that the location of the first micro-crack is sub-surface which is consistent with observations made by a number of

Fig. 9. (a) Discrete representation of the semi-infinite domain simulating a rolling contact; (b) zoomed view.

354

N. Raje et al. / International Journal of Fatigue 31 (2009) 346–360

researchers for lubricated rolling contacts with smooth surfaces [44–47]. Note that the friction coefficient used in this study is representative of these lubrication and surface conditions. The initiation life is defined as the number of cycles elapsed before the first micro-crack appears. Subsequently, multiple joints are damaged forming several micro-cracks. The micro-cracks are primarily found to be oriented parallel to the surface which is consistent with experimental observations [48,49]. With further cycling, some of the micro-cracks coalesce into longer cracks and form damaged zones with distinct trajectories. The damaged zones then propagate towards the contact surface and form an elliptical spall. The subsequent number of cycles after crack initiation before a surface crack forms constitute the propagation life. The computation is stopped once a surface breaking crack forms. This constitutes the total life. It is to be recognized that the entire process is complex and involves considerable inter-play between developing sets of micro-cracks. Experimental studies on the continuous observation of subsurface damage in rolling contacts have been few because of the difficulty in sectioning the specimen at the right location. Chen et al. [46] carried out a study to observe the continuous evolution of sub-surface crack patterns in lubricated rolling line contacts for Chinese bearing steel GCr15 which is similar in composition and hardness to bearing steel AISI-52100 used in the present analysis. Their results are shown in Fig. 10b. The results obtained from our model are found to be in good qualitative agreement with these experimental observations. Note that the long crack that grows into the material in the experimental images is due to the effect of circumferential residual tensile stresses that develop below the bearing race. Such an effect is absent in the present analysis since plasticity is not incorporated in the model. Hence, this crack is not observed in the simulation results. Fig. 11 shows damage evolution in a randomly selected interelement joint in the sample domain for damage increments of DD = 0.1 and 0.2. The damage evolution curves for the two damage increments are close to one another however, the computation time for the lower damage increment is twice that for the higher increment. Fig. 12 shows the relative proportions of lives spent in the initiation and propagation stages. Note that in this case, a significant fraction of the total life is spent in the crack propagation stage. This is primarily due to the spatially varying shear stress field existing in rolling contacts. Shear stress is maximum at some depth below the surface and decreases in the direction of the surface. Cracks are generated in the high shear stress region subsurface. As they progress towards the surface they encounter diminishing magnitudes of shear stresses, which leads to retardation in their growth rates. As a result, it takes a considerable amount of time before a sub-surface crack reaches the surface. There is substantial experimental evidence in literature that subsurface cracks are formed early in the fatigue life of rolling contacts and the major proportion of life is spent in the crack propagation stage [46,47]. The results from the current model are in good qualitative agreement with experimental observations. 7.2. Effect of topological material disorder To study only the effect of material topology (microstructural geometric randomness) on the fatigue life, material properties (E, rr and m) are kept uniform throughout the domain. Analyses are carried out for the 40 material domains using a maximum Hertzian pressure of pmax = 2 GPa and damage increments of DD = 0.1 and 0.2. Fig. 13 depicts the Weibull life plots for the two cases. The two damage increments yield initiation lives that are about 19% different, however the total lives are within 5% of each other, with the smaller increment yielding lower lives. However, considering the computational efficiency achievable by using a larger dam-

Fig. 10a. Evolution of sub-surface damage in rolling contacts through micro-crack initiation and coalescence.

age increment, all subsequent analyses are performed using DD = 0.2. The initiation lives do not show much scatter between domains and the overall scatter in the total lives is primarily governed by the scatter in propagation lives. A number of statistical distributions can be potentially used to fit the simulation data, e.g. normal,

N. Raje et al. / International Journal of Fatigue 31 (2009) 346–360

e  

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

355

ð38Þ

Here, F is the probability that life is less than some number of cycles N, while a P 0, b > 0 and e > 0 are the three parameters for the distribution. Note that a is the minimum life parameter. In the absence of this parameter (a = 0), the life distribution takes the form of the popular 2-parameter Weibull distribution and e becomes the slope for the linear plot (Weibull slope). The Weibull parameters are obtained using the method of Maximum Likelihood Estimate (MLE) [51] and are listed in Table 2. The goodness-of-fit for the two distributions is assessed according to the Anderson-Darling test (AD value) and the corresponding P value. Lower the AD value and higher the P value, better is the fit. It is seen from Table 2 that the 3parameter Weibull distribution gives a much better fit compared to the 2-parameter distribution. A life quantity that is commonly used in rolling contacts is the life at 90% reliability or 10% probability of failure also known as the L10 life. This life quantity is used for demonstrating the results of the present analysis. 7.3. Effect of variable elastic modulus and initial flaws

Fig. 10b. Evolution of sub-surface crack patterns in rolling contacts as observed under microscope [46].

lognormal, Weibull, gamma, exponential. However, lives of rolling element bearings that experience rolling contact stressing are found to follow the Weibull distribution closely [41]. Hence, this distribution is used to fit the simulation data. Weibull slopes (e) for the initiation and total lives using a 2-parameter Weibull fit are 11.75 and 1.85, respectively. Experimentally obtained Weibull slopes for rolling element bearings range from 1.0–3.5 [50]. Also, the expected Weibull slope for rolling element bearing lives as per several manufacturers’ catalogs is 1.5. Note that these slopes correspond to the total failure lives. It needs to be pointed out that results for the total lives show a small curvature, especially at low probabilities of failure. This alludes to the fact that a minimum life exists even for 100% reliability which physically relates to the minimum amount of time required for crack(s) to generate in the sub-surface region, propagate to the surface and form a spall. In such circumstances, a 3-parameter Weibull fit with a minimum life parameter is a better approximation for the resulting life distribution. The 3-parameter Weibull distribution is given by,

A non-uniform distribution of elastic modulus can occur because of the presence of different material phases or inclusions within the base material matrix. This results in local notch effects due to strain mismatch and assists the fatigue failure process. In this investigation, the elastic modulus E is assumed to vary spatially, following a Normal distribution centered around 200 GPa such that: 180 GPa < E < 220 GPa. Fig. 14 illustrates the distribution of lives obtained for the 40 material domains. The total life distributions again follow a 3-parameter characteristic due to the curved nature of the results. Table 2 depicts results obtained using statistical analyses. From the results, it can be inferred that there is no significant effect of the non-uniform elastic modulus on the overall spalling lives. However, scatter in lives increases (Weibull slope e decreases) with increase in the range of elastic modulus variation. Flaws present in the material act as likely sites for crack initiation. Here, the effect of 2 initial flaws on fatigue lives is studied. Flaws are modeled as failed joints (i.e., joints with damage variable D = 1). The elastic modulus is kept uniform (E = 200 GPa) throughout the domain. The flaw location and orientation is randomly selected for every material domain. The resultant life distributions for the 40 domains are depicted in Fig. 14 and the results are listed in Table 2. It can be seen that the flaws do not affect the total lives significantly; however, they do affect the initiation lives for the next micro-crack to form. Overall, the scatter in lives increases and the L10 life decreases slightly with increasing number of flaws. Again, from the AD and P values, it is seen that the 3-parameter Weibull distribution fits the data better. 7.4. Effect of stress level In order to study the effect of stress level on fatigue lives in rolling contacts, the Hertzian pressure pmax is varied from 1 GPa to 2.5 GPa. Note that this is the nominal range of contact pressures experienced in line contact roller bearings. Table 3 lists the Weibull parameters obtained from statistical analyses. The 3-parameter Weibull distribution is found to fit the data better than the 2-parameter Weibull distribution for the entire range of pressures considered. Fig. 15 shows the variation of L10 life with contact pressure. Experimental endurance data for cylindrical roller bearings made from carbon vacuum-degassed AISI-52100 bearing steel [52] has been included for comparison. It is to be noted that contacts in cylindrical roller bearings are close to the ideal line contact scenario assumed here and offer a good comparison between the

356

N. Raje et al. / International Journal of Fatigue 31 (2009) 346–360

Fig. 11. Damage evolution for a randomly selected inter-element joint under rolling contact loading.

Fig. 12. Progression of stress cycles with loading blocks for rolling contact.

present model and experiments. The model is found to fit within the range of experimental data. Note that there is considerable scatter in the experimental data since it has been compiled from disparate sources and under different surface lubrication conditions (e.g. mixed or boundary lubrication). In the present analysis, the surface pressure is assumed to be perfectly Hertzian and surfaces are assumed to be perfectly smooth. A more justifiable comparison can be made with the Lundberg-Palmgren theory [41] which forms the basis of current life ratings for rolling element bearings. The comparison has been depicted in Fig. 15 for a 209

series cylindrical roller bearing having line contacts. Good corroboration is found between the model and the Lundberg-Palmgren theory. The degree of scatter in fatigue lives is found to be fairly consistent over the range of contact pressures considered. 7.5. Relative initiation and propagation lives Fig. 16 shows the percentage of life spent in the crack initiation stage as a function of total life and contact pressure. The initiation and propagation lives are the average for the 40 material domains.

357

N. Raje et al. / International Journal of Fatigue 31 (2009) 346–360

Fig. 13. Weibull life plots for different damage increments DD.

Table 2 Weibull statistics for the 40 material domains under rolling contact loading (pmax = 2 GPa) Material condition

Baseline 180 < E < 220 GPa 2 flaws

2-parameter Weibull

3-parameter Weibull

e

L10 ( 106)

AD

P

e

a ( 106)

b ( 106)

L10 ( 106)

AD

P

1.85 1.61 1.81

23.47 20.75 22.71

0.80 0.68 0.83

0.04 0.07 0.03

1.27 1.17 1.31

18.11 14.99 15.78

55.99 63.11 58.07

27.64 24.15 26.22

0.38 0.33 0.42

0.42 >0.5 0.35

Fig. 14. Weibull life plots for varying material conditions.

358

N. Raje et al. / International Journal of Fatigue 31 (2009) 346–360

Table 3 Weibull statistics for varying contact pressures (E = 200 GPa) pmax (GPa)

1.0 1.3 1.5 1.8 2.0 2.2 2.5

2-parameter Weibull

3-parameter Weibull 6

e

L10 ( 10 )

AD

P

e

a ( 106)

b ( 106)

L10 ( 106)

AD

P

1.97 1.82 1.96 2.01 1.85 1.73 1.52

15604 1168.4 313.36 61.07 23.47 10.19 3.32

1.42 1.45 1.47 1.13 0.80 0.66 0.57

<0.01 <0.01 <0.01 <0.01 0.04 0.08 0.14

1.22 1.15 1.27 1.38 1.27 1.22 1.18

15823 1214.7 288.12 47.49 18.11 7.49 2.06

29284 2453.7 621.84 129.02 55.99 27.35 11.65

20478 1561.9 393.34 72.83 27.64 11.85 3.78

0.45 0.45 0.58 0.56 0.38 0.34 0.35

0.3 0.29 0.14 0.15 0.42 >0.5 0.49

Fig. 15. Variation of L10 life with rolling contact pressure.

Fig. 16. Percentage of life spent in the crack initiation stage.

N. Raje et al. / International Journal of Fatigue 31 (2009) 346–360 Table 4 Comparison of model with experimental data for percentage of life spent in crack initiation stage Source

pmax (GPa)

% initiation life

Model Chen et al. [46]

2.5 2.45

12.05 10

As expected, at higher contact pressures (low total life), cracks are formed relatively early and the crack initiation stage is a small fraction of the total life. As the contact pressure is decreased, the percentage of life spent in the initiation stage gradually increases. The model results are compared with available experimental data from Chen et al. [46] who studied the relative proportions of fatigue life spent in initiation and propagation phases for rolling line contacts. A comparison between the experimental results of Chen et al. [46] and the current model is given in Table 4. It is interesting to note that in the experimental study by Chen et al. [46], the shortest visible crack length defining the crack initiation stage was 10–15 lm which is roughly the size of a micro-crack considered in the present analysis and the length of a grain. From Table 4, it is seen that model results are in good quantitative agreement with experiments. 8. Summary and conclusions This paper presents a stochastic fatigue model to study the effect of microstructural topological disorder on fatigue damage initiation and propagation in polycrystalline materials. The model incorporates the principles of damage mechanics in a discrete material framework and is capable of simulating the inherent microstructural disorder that leads to fatigue life scatter. Three distinct sources are found to cause dispersion in fatigue lives: (i) the topological (geometric) variation in the material microstructure, (ii) inhomogeneities due to variation of material properties throughout the microstructure and (iii) presence of internal flaws. Micro-crack initiation, coalescence and propagation stages are modeled using damaged zones in a unified framework. The model is applied to study fatigue in rolling contacts. Fatigue lives obtained from the model are in good agreement with experimental observations. The propagation life constitutes a significant fraction of the total life. This is consistent with experimental observations. Moreover, the lives follow a 3-parameter Weibull distribution. The relative proportion of lives spent in the initiation and propagation stages are in good quantitative agreement with experiments. Acknowledgement The authors would like to express their deepest appreciation to the Honeywell Company for their support of this project. References [1] Schlangen E, Garboczi EJ. New method for simulating fracture using an elastically uniform random geometry lattice. Int J Eng Sci 1996;34(10): 1131–44. [2] Rinaldi A, Peralta P, Krajcinovic D, Lai YC. Prediction of scatter in fatigue properties using discrete damage mechanics. Int J Fatigue 2006;28:1069–80. [3] Xu XP, Needleman A. Numerical simulation of dynamic interfacial crack growth allowing for crack growth away from the bond line. Int. J. Fract 1995;74:253–75. [4] Camacho G, Ortiz M. Computational modeling of impact damage in brittle materials. Int J Solid Struct 1996;33:2899–938. [5] Espinosa H, Zavattieri P, Dwivedi S. A finite deformation continuum/discrete model for the description of fragmentation and damage in brittle materials. J Mech Phys Solid 1998;46(10):1909–42.

359

[6] Miller O, Freund L, Needleman A. Modeling and simulation of dynamic fragmentation in brittle materials. Int J Fract 1999;96(2):101–25. [7] Bolander JE, Saito S. Fracture analysis using spring networks with random geometry. Eng Fract Mech 1998;61(5-6):569–91. [8] Kawai T. New discrete models and their application to seismic response analysis of structures. Nucl Eng Design 1978;48:207–29. [9] Herrmann HJ, Roux S. Statistical Models For The Fracture of Disordered Media. Amsterdam: Elsevier/North Holland; 1990. [10] Curtin WA, Scher H. Brittle fracture in disordered materials: a spring network model. J Mater Res 1990;5(3):535–53. [11] Jagota A, Bennison SJ. Spring-network and finite element models for elasticity and fracture. In: Bardhan KK, Chakrabarti BK, Hansen A, editors. Non-Linearity and Breakdown in Soft Condensed Matter (Springer Lecture Notes in Physics 437). Berlin: Springer-Verlag; 1994. p. 186–201. [12] Espinosa HD, Zavattieri PD. A grain level model for the study of failure initiation and evolution in polycrystalline brittle materials. Part I: theory and numerical implementation. Mech Mater 2003;35:333–64. [13] Herrmann HJ, Hansen A, Roux S. Fracture of disordered elastic lattices in two dimensions. Phys Rev B 1989;39(1):637–48. [14] Meakin GL, Sander LM, Louis E, Guinea F. A simple two-dimensional model for crack propagation. J Phys A: Math Gen 1989;22:1393–403. [15] Bazant ZP, Tabbara MR, Kazemi MT, Pijaudier-Cabot G. Random particle model for fracture of aggregate or fiber composites. J Eng Mech 1990;116(8): 1686–705. [16] de Andres A, Perez J, Ortiz M. Elastoplastic finite element analysis of threedimensional fatigue crack growth in aluminum shafts subjected to axial loading. Int J Solid Struct 1999;36:2231–58. [17] Nguyen O, Repetto EA, Ortiz M, Radovitzky RA. A cohesive model of fatigue crack growth. Int J Fract 2001;110:351–69. [18] Roe KL, Siegmund T. An irreversible cohesive zone model for interface fatigue crack growth simulation. Eng Fract Mech 2003;70:209–32. [19] Yang B, Mall S, Ravi-Chander K. A cohesive zone model for fatigue crack growth in quasi brittle materials. Int J Solid Struct 2001;38:3927–44. [20] Iesulauro E, Ingraffea AR, Arwade S, Wawrzynek PA. Simulation of grain boundary decohesion and crack initiation in aluminum microstructure models. ASTM STP 2002;1417:715–28. [21] Lemaitre J. A Course on Damage Mechanics. Berlin: Springer-Verlag; 1992. [22] Raje NN, Sadeghi F, Rateick Jr RG. A discrete element approach to evaluate stresses due to line loading on an elastic half-space. Comp Mech 2007;40(3):513–29. [23] Potapov AV, Campbell CS, Hopkins MA. A two-dimensional dynamic simulation of solid fracture. Part I: description of the model. Int J Mod Phys C 1995;6(3):371–98. [24] Andersson J. The influence of grain size variation on metal fatigue. Int J Fatigue 2005;27:847–52. [25] Ito O, Fuller ER. Computer modeling of anisotropic grain microstructure in two dimensions. Acta Metall Mater 1993;41(1):191–8. [26] Zavattieri PD, Espinosa HD. Grain level analysis of crack initiation and propagation in brittle materials. Acta Mater 2001;49:4291–311. [27] Preparata FP, Shamos MI. Computational Geometry: an introduction. Berlin: Springer-Verlag; 1985. [28] L.M. Kachanov Time of the rupture process under creep conditions. Izv Akad Nauk SSR, Otd Tekh Nauk 8 (1958) 26-31. [29] Robotnov YN. Creep Problems in Structural Mechanics. Amsterdam: NorthHolland; 1969. [30] Bolotin VV, Belousov IL. Early fatigue crack growth as the damage accumulation process. Probab Eng Mech 2001;16:279–87. [31] Chaboche JL, Lesne PM. A non-linear continuous fatigue damage model. Fatigue Fract Eng Mater Struct 1988;11(1):1–17. [32] Xiao YC, Li S, Gao Z. A continuum damage mechanics model for high cycle fatigue. Int J Fatigue 1998;20(7):503–8. [33] Memon IR, Cui D, Zhang X. Fatigue life prediction of 3-D problems by damage mechanics – finite element additional load method. In: Wu XR, Wang ZG, editors. Fatigue ‘99: Proceedings of the Seventh International Fatigue Congress, vol. 2. Beijing: Higher Education Press; 1999. p. 827–32. [34] Bolotin VV. Mechanics of Fatigue. Boca Raton: CRC Press; 1999. [35] Raje N, Sadeghi F, Rateick Jr RG, Hoeprich MR. A numerical model for life scatter in rolling element bearings. J Trib 2008;130(1):0110011. [36] Lemaitre J, Sermage JP, Desmorat R. A two scale damage concept applied to fatigue. Int J Fract 1998;97:67–81. [37] Cundall PA, Strack ODL. A discrete numerical model for granular assemblages. Geotech 1979;29:47–65. [38] Hocking G. The discrete element method for analysis of fragmentation of discontinua. Eng Comp 1992;9:145–55. [39] M.A. Hopkins The numerical simulation of systems of multitudinous polygonal blocks. US Army Cold Regions Research and Engineering Laboratory, USACRREL Report CR 99-22; 1992. [40] Styri H. Fatigue strength of ball bearing races and heat-treated 52100 steel specimens. Proc ASTM 1951;51:682–700. [41] Lundberg G, Palmgren A. Dynamic capacity of rolling bearings. Acta Polytech Mech Eng Ser 1947;1(3):7–50. [42] W.E. Littmann The mechanism of contact fatigue. Interdisciplinary approach to the lubrication of concentrated contacts. Ku PM, editor. Proceedings of a Symposium, Troy, NY, NASA Special Report, SP-237 1969, pp. 309-378. [43] Johnson KL. Contact Mechanics. Cambridge: Cambridge University press; 1985.

360

N. Raje et al. / International Journal of Fatigue 31 (2009) 346–360

[44] Shao E, Huang X, Wang C, Zhu Y, Chen Q. A method of detecting rolling contact crack initiation and the establishment of crack propagation curves. STLE Transactions 1987;31(1):6–11. [45] Leng X, Chen Q, Shao E. Initiation and propagation of case crushing cracks in rolling contact fatigue. Wear 1988;122:33–43. [46] 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–18. [47] Nelias D, Dumont ML, Champiot F, Vincent A, Girodin D, Fougeres R, et al. Role of inclusions, surface roughness and operating conditions on rolling contact fatigue. J Trib 1999;121(2):240–51.

[48] Yoshioka T. Detection of rolling contact sub-surface fatigue cracks using acoustic emission technique. Lub Eng 1993;49(4):303–8. [49] Otsuka A, Sugawara H, Shomura M. A test method for mode II fatigue crack growth relating to a model for rolling contact fatigue. Fatigue Fract Eng Mater Struct 1996;19(10):1265–75. [50] Harris TA. Rolling Bearing Analysis. New York: John Wiley & Sons; 2001. [51] Lockhart RA, Stephens MA. Estimation and tests of fit for the threeparameter Weibull distribution. J Royal Stat Soc B 1994;56(3):491–500. [52] Harris TA, Barnsby RM. Life ratings for ball and roller bearings. Proc IMechE J: J Eng Trib 2001;215(6):577–95.