Fracture Mechanics Based Models of Structural and Contact Fatigue

The subsurface initiated contact fatigue failure is one of the dominating mechanisms of failure of moving machine parts involved in cyclic motion. Structural fatigue failure may be of surface or subsurface origin. The analysis of a significant amount of accumulated experimental data obtained from field exploitation and laboratory testing provides undisputable evidence of the most important factors affecting contact fatigue [1]. It is clear that the factors affecting contact fatigue the most are as follows (a) acting cyclic normal stress and frictional stress (detailed lubrication conditions, surface roughness, etc.) which in part are determined by the part geometry, (b) distribution of residual stress versus depth, (c) initial statistical defect/crack distribution versus defect size, and location, (d) material elastic and fatigue parameters as functions of materials hardness, etc., (e) material fracture toughness, (f) material hardness versus depth, (g) machining and finishing operations, (h) abrasive contamination of lubricant and residual surface contamination, (i) non-steady cyclic loading regimes, etc. In case of structural fatigue the list of the most important parameters affecting fatigue performance is similar. None of the existing contact or structural fatigue models developed for prediction of contact fatigue life of bearings and gears as well as other structures takes into account all of the above operational and material conditions. Moreover, at best, most of the existing contact fatiguemodels are only partially based on the fundamental physical and mechanical mechanisms governing the fatigue phenomenon. Most of these models are of empirical nature and are based on assumptions some ofwhich are not supported by experimental data or are controversial as it is in the case of fatigue models for bearings and gears [1]. Some models involve a number of approximations that usually do not reflect the actual processes occurring in material.


Introduction
The subsurface initiated contact fatigue failure is one of the dominating mechanisms of failure of moving machine parts involved in cyclic motion. Structural fatigue failure may be of surface or subsurface origin.
The analysis of a significant amount of accumulated experimental data obtained from field exploitation and laboratory testing provides undisputable evidence of the most important factors affecting contact fatigue [1]. It is clear that the factors affecting contact fatigue the most are as follows (a) acting cyclic normal stress and frictional stress (detailed lubrication conditions, surface roughness, etc.) which in part are determined by the part geometry, (b) distribution of residual stress versus depth, (c) initial statistical defect/crack distribution versus defect size, and location, (d) material elastic and fatigue parameters as functions of materials hardness, etc., (e) material fracture toughness, (f) material hardness versus depth, (g) machining and finishing operations, (h) abrasive contamination of lubricant and residual surface contamination, (i) non-steady cyclic loading regimes, etc. In case of structural fatigue the list of the most important parameters affecting fatigue performance is similar. None of the existing contact or structural fatigue models developed for prediction of contact fatigue life of bearings and gears as well as other structures takes into account all of the above operational and material conditions. Moreover, at best, most of the existing contact fatigue models are only partially based on the fundamental physical and mechanical mechanisms governing the fatigue phenomenon. Most of these models are of empirical nature and are based on assumptions some of which are not supported by experimental data or are controversial as it is in the case of fatigue models for bearings and gears [1]. Some models involve a number of approximations that usually do not reflect the actual processes occurring in material.
Therefore, a comprehensive mathematical models of contact and structural fatigue failure should be based on clearly stated mechanical principles following from the theory of elasticity, lubrication theory of elastohydrodynamic contact interactions, and fracture mechanics. Such models should take into consideration all the parameters described in items (a)-(e) and beyond. The advantage of such comprehensive models would be that the effect of variables such as steel cleanliness, externally applied stresses, residual stresses, etc. on contact and structural fatigue life could be examined as single or composite entities.
The goal of this chapter is to provided fracture mechanics based models of contact and structural fatigue. Historically, one of the most significant problems in realization of such an approach is the availability of simple but sufficiently precise solutions for the crack stress intensity factors. To overcome this difficulty some problems of fracture mechanics will be analyzed and their solutions will be represented in an analytical form acceptable for the further usage in modeling of contact and structural fatigue. In particular, a problem for an elastic half-plain weakened with a number of subsurface cracks and loaded with contact normal and frictional stresses as well with residual stress will be formulated and its asymptotic analytical solution will be presented. The latter solutions for the crack stress intensity factors are expressed in terms of certain integrals of known functions. These solutions for the stress intensity factors will be used in formulation of a two-dimensional contact fatigue model. In addition to that, a three-dimensional model applicable to both structural and contact fatigue will be formulated and some examples will be given. The above models take into account the parameters indicated in items (a)-(e) and are open for inclusion of the other parameters significant for fatigue.
In particular, these fatigue models take into account the statistical distribution of inclusions/cracks over the volume of the material versus their size and the resultant stress acting at the location of every inclusion/crack. Some of the main assumptions of the models are that (1) fatigue process in any machine part or structural unit runs in a similar manner and it is a direct reflection of the acting cyclic stresses and material properties and (2) the main part of fatigue life corresponds to the fatigue crack growth period, i.e. fatigue crack initiation period can be neglected. Furthermore, the variation in the distribution of cracks over time due to their fatigue growth is accounted for which is absent in all other existing fatigue models. The result of the above fatigue modeling is a simple relationship between fatigue life and cyclic loading, material mechanical parameters and its cleanliness as well as part geometry. are observed. In such cases these defect clusters can be represented by single defects of approximately the same size. Suppose there is a characteristic size L σ in material that is determined by the typical variations of the material stresses, grain and surface geometry. It is also assumed that there is a size L f in material such that L d L f L σ , where L d is the typical distance between the material defects. In other words, it is assumed that the defect population in any such volume L 3 f is large enough to ensure an adequate statistical representation. It means that any parameter variations on the scale of L f are indistinguishable for the fatigue analysis purposes and that in the further analysis any volume L 3 f can be represented by its center point (x, y, z).
Therefore, there is an initial statistical defect distribution in the material such that each defect can be replaced by a subsurface penny-shaped or a surface semi-circular crack with a radius approximately equal to the half of the defect diameter. The usage of penny-shaped subsurface and semi-circular surface cracks is advantageous to the analysis because such fatigue cracks maintain their shape and their size is characterized by just one parameter. The orientation of these crack propagation will be considered later. The initial statistical distribution is described by a probabilistic density function f (0, x, y, z, l 0 ), such that f (0, x, y, z, l 0 )dl 0 dxdydz is the number of defects with the radii between l 0 and l 0 + dl 0 in the material volume dxdydz centered about point (x, y, z). The material defect distribution is a local characteristic of material defectiveness. The model can be developed for any specific initial distribution f (0, x, y, z, l 0 ). Some experimental data [3] suggest a log-normal initial defect distribution f (0, x, y, z, l 0 ) versus the defect initial radius l 0 where μ ln and σ ln are the mean value and standard deviation of the crack radii, respectively.

Direction of fatigue crack propagation
It is assumed that the duration of the crack initiation period is negligibly small in comparison with the duration of the crack propagation period. It is also assumed that linear elastic fracture mechanics is applicable to small fatigue cracks. The details of the substantiation of these assumptions can be found in [1]. Based on these assumptions in the vicinity of a crack the stress intensity factors completely characterize the material stress state. The normal k 1 and shear k 2 and k 3 stress intensity factors at the edge of a single crack of radius l can be represented in the form [4] where σ 1 is the maximum of the local tensile principal stress, F 1 , F 2 , and F 3 are certain functions of the point coordinates (x, y, z) and the crack orientation angles α and β with respect to the coordinate planes. The coordinate system is introduced in such a way that the x-and y-axes are directed along the material surface while the z-axis is directed perpendicular to the material surface.
The resultant stress field in an elastic material is formed by stresses σ x (x, y, z), σ y (x, y, z), σ z (x, y, z), τ xz (x, y, z), τ xy (x, y, z), and τ zy (x, y, z). Some regions of material are subjected to tensile stress while other regions are subjected to compressive stress. Conceptually, there is no difference between the phenomena of structural and contact fatigue as the local material response to the same stress in both cases is the same and these cases differ in their stress fields only. As long as the stress levels do not exceed the limits of applicability of the quasi-brittle linear fracture mechanics when plastic zones at crack edges are small the rest of the material behaves like an elastic solid. The actual stress distributions in cases of structural and contact fatigue are taken in the proper account. In contact interactions where compressive stress is usually dominant there are still zones in material subjected to tensile stress caused by contact frictional and/or tensile residual stress [1,5].
Experimental and theoretical studies show [1] that after initiation fatigue cracks propagate in the direction determined by the local stress field, namely, perpendicular to the local maximum tensile principle stress. Therefore, it is assumed that fatigue is caused by propagation of penny-shaped subsurface or semi-circular surface cracks under the action of principal maximum tensile stresses. Only high cycle fatigue phenomenon is considered here. On a plane perpendicular to a principal stress the shear stresses are equal to zero, i.e. the shear stress intensity factors k 2 = k 3 = 0. To find the plane of fatigue crack propagation (i.e. the orientation angles α and β), which is perpendicular to the maximum principal tensile stress, it is necessary to find the directions of these principal stresses. The latter is equivalent to solving the equations k 2 (N, α, β, l, x, y, z) = 0, k 3 (N, α, β, l, x, y, z) = 0.
Usually, there are more than one solution sets to these equations at any point (x, y, z). To get the right angles α and β one has to chose the solution set that corresponds to the maximum tensile principal stress, i.e. maximum of the normal stress intensity factor k 1 (N, l, x, y, z). That guarantees that fatigue cracks propagate in the direction perpendicular to the maximum tensile principal stress if α and β are chosen that way.
For steady cyclic loading for small cracks k 20 = k 2 / √ l and k 30 = k 3 / √ l are independent from the number of cycles N and crack radius l together with equation (3) lead to the conclusion that for cyclic loading with constant amplitude the angles α and β characterizing the plane of fatigue crack growth are independent from N and l. Thus, angles α and β are functions of only crack location, i.e. α = α(x, y, z) and β = β(x, y, z). For the most part of their lives fatigue cracks created and/or existed near material defects remain small. Therefore, penny-shaped subsurface cracks conserve their shape but increase in size.
Even for the case of an elastic half-space it is a very difficult task to come up with sufficiently precise analytical solutions for the stress intensity factors at the edges of penny-shaped subsurface or semi-circular surface crack of arbitrary orientation at an arbitrary location (x, y, z). However, due to the fact that practically all the time fatigue cracks remain small and exsert little influence on the material general stress state the angles α and β can be determined in the process of calculation of the maximum tensile principle stress σ. The latter is equivalent to solution of equations (3). As soon as σ is determined for a subsurface crack its normal stress intensity factor k 1 can be approximated by the normal stress intensity factor for the case of a single crack of radius l in an infinite space subjected to the uniform tensile stress σ, i.e. by k 1 = 2σ √ l/π.

Fatigue crack propagation
Fatigue cracks propagate at every point of the material stressed volume V at which max where the maximum is taken over the duration of the loading cycle T and k th is the material stress intensity threshold. There are three distinct stages of crack development: (a) growth of small cracks, (b) propagation of well-developed cracks, and (c) explosive and, usually, unstable growth of large cracks. The stage of small crack growth is the slowest one and it represents the main part of the entire crack propagation period. This situation usually causes confusion about the duration of the stages of crack initiation and propagation of small cracks. The next stage, propagation of well-developed cracks, usually takes significantly less time than the stage of small crack growth. And, finally, the explosive crack growth takes almost no time.
A relatively large number of fatigue crack propagation equations are collected and analyzed in [6]. Any one of these equations can be used in the model to describe propagation of fatigue cracks. However, the simplest of them which allows to take into account the residual stress and, at the same time, to avoid the usage of such an unstable characteristic as the stress intensity threshold k th is Paris's equation where g 0 and n are the parameters of material fatigue resistance and l 0 is the crack initial radius. Notice, that in cases of loading and relaxation such as in contact fatigue k 1 = k 1 .
Fatigue cracks propagate until they reach their critical size with radius l c for which k 1 = K f (K f is the material fracture toughness, i.e. to the radius of l c = (K f /k 10 ) 2 ). After that their growth becomes unstable and very fast. Usually, the stage of explosive crack growth takes just few loading cycles.
It can be shown that the number of loading cycles needed for a crack to reach its critical radius is almost independent from the material fracture toughness K f . This conclusion is supported by direct numerical simulations. For the further analysis, it is necessary to determine for a crack its initial radius l 0c , which after N loading cycles reaches the critical size of l c . Solving the initial-value problem (4) one obtains the formula where l 0c depends on N, x, y, and z. Obviously, for n > 2 and fixed x, y, and z the value of l 0c is a decreasing function of N. It is important to keep in mind that l 0c (N, x, y, z) is minimal where k 10 (x, y, z) is maximal, which, in turn, happens where the material tensile stress reaches its maximum.

Crack propagation statistics
To describe crack statistics after the crack initiation stage is over it is necessary to make certain assumptions. The simplest assumptions of this kind are: the existing cracks do not heal and new cracks are not created. In other words, the number of cracks in any material volume remains constant in time. Based on a practically correct assumption that the defect distribution is initially scarce, the coalescence of cracks and changes in the general stress field are possible only when cracks have already reached relatively large sizes. However, this may happen only during the last stage of crack growth the duration of which is insignificant for calculation of fatigue life. Therefore, it can be assumed that over almost all life span of fatigue cracks their orientations do not change. This leads to the equation for the density of crack distribution f (N, x, y, z, l) as a function of crack radius l after N loading cycles in a small parallelepiped dxdydz with the center at the point with coordinates (x, y, z) which being solved for f (N, x, y, z, l) gives where l 0 and dl 0 /dl as functions of N and l can be obtain from the solution of (4) in the form Equations (7) and (8) lead to the expression for the crack distribution function f after N loading cycles f (N, x, y, z, l) = f (0, x, y, z, l 0 (N, l, y, z)){1 where l 0 (N, l, y, z) is determined by the first of the equations in (8).
Formula (9) leads to a number of important conclusions. The fatigue crack distribution function f (N, x, y, z, l) depends on the initial crack distribution f (0, x, y, z, l 0 ) and it changes with the number of applied loading cycles N in such a way that the crack volume density ρ(N, x, y, z) remains constant. Because of crack growth the crack distribution f (N, x, y, z, l) widens with respect to l with number of loading cycles N.

Local fatigue damage accumulation
Let us design the measure of material fatigue damage. If at a certain point (x, y, z) after N loading cycles radii of all cracks l < l c then there is no damage at this point and the material local survival probability p(N, x, y, z) = 1. On the other hand, if at this point after N loading cycles radii of all cracks l ≥ l c , then all cracks reached the critical size and the material at this point is completely damaged and the local survival probability p(N, x, y, z) = 0. Obviously, the more fatigue cracks with larger radii l exist at the point the lower is the local survival probability p(N, x, y, z). It is reasonable to assume that the material local survival probability p(N, x, y, z) is a certain monotonic measure of the portion of cracks with radius l below the critical radius l c . Therefore, p(N, x, y, z) can be represented by the expressions Obviously, the local survival probability p(N, x, y, z) is a monotonically decreasing function of the number of loading cycles N because fatigue crack radii l tend to grow with the number of loading cycles N.
To calculate p(N, x, y, z) from (10) one can use the specific expression for f determined by (9). However, it is more convenient to modify it as follows where l 0c is determined by (5) and ρ is the initial volume density of cracks. Thus, to every material point (x, y, z) is assigned a certain local survival probability p(N, x, y, z), 0 ≤ p(N, x, y, z) ≤ 1.
Equations (11) demonstrate that the material local survival probability p(N, x, y, z) is mainly controlled by the initial crack distribution f (0, x, y, z, l 0 ), material fatigue resistance parameters g 0 and n, and external contact and residual stresses. Moreover, the material local survival probability p(N, x, y, z) is a decreasing function of N because l 0c from (5) is a decreasing function of N for n > 2.

Global fatigue damage accumulation
The survival probability P(N) of the material as a whole is determined by the local probabilities of all points of the material at which fatigue cracks are present. It is assumed that the material fails as soon as it fails at just one point. It is assumed that the initial crack distribution in the material is discrete.
where N c is the total number of points in the material stressed volume V at which fatigue cracks are present. Then based on the above assumption the material survival probability P(N) is equal to Obviously, probability P(N) from (12) satisfies the inequalities In (13) the right inequality shows that the survival probability P(N) is never greater than the minimum value p m (N) of the local survival probability p(N, x, y, z) over the material stressed volume V.
An analytical substantiation for the assumption that the first pit is created by the cracks from a small material volume with the smallest survival probability p m (N) is provided in [1]. Moreover, the indicated analysis also validates one of the main assumptions of the model that new cracks are not being created. Namely, if new cracks do get created in the process of loading, they are very small and have no chance to catch up with already existing and propagating larger cracks. The graphical representation of this fact is given in Fig. 1 [1]. In this figure fatigue cracks are initially randomly distributed over the material volume with respect to their normal stress intensity factor k 1 and are allowed to grow according to Paris' law (see (4)) with sufficiently high value of n = 6.67 − 9. In Fig. 1 the values of the normal stress intensity factor k 1 are shown at different time moments (k 0 and L 0 are the characteristic normal stress intensity factor and geometric size of the solid). These graphs clearly show that a crack with the initially larger value of the normal stress intensity factor k 1 propagates much faster than all other cracks, i.e. the value of its k 1 increases much faster than the values of k 1 for all other cracks, which are almost dormant. As a result of that, the crack with the initially larger value of k 1 reaches its critical size way ahead of other cracks. This event determines the time and the place where fatigue failure occurs initially. Therefore, in spite of formula (12) which indicates that all fatigue cracks have influence on the survival probability P(N), for high values of n the material survival probability P(N) is a local fatigue characteristic, and it is determined by the material defect with the initially highest value of the stress intensity factor k 1 . The higher the power n is the more accurate this approximation is.
Therefore, assuming that it is a very rare occurrence when more than one fatigue initiation/spall happen simultaneously, it can be shown that at the early stages of the fatigue process the material global survival probability P(N) is determined by the minimum of the local survival probability p m (N) where the maximum is taken over the (stressed) volume V of the solid.
If the initial crack distribution is taken in the log-normal form (1) then where er f (x) is the error integral [8]. Obviously, the local survival probability p m (N) is a complex combined measure of applied stresses, initial crack distribution, material fatigue parameters, and the number of loading cycles.
In cases when the mean μ ln and the standard deviation σ ln are constants throughout the material formula (15) can be significantly simplified To determine fatigue life N of a contact for the given survival probability P(N) = P * , it is necessary to solve the equation

Fatigue life calculation
Suppose the material failure occurs at point (x, y, z) with the probability 1 − P(N). That actually determines the point where in (16) the minimum over the material volume V is reached. Therefore, at this point in (16), the operation of minimum over the material volume V can be dropped. By solving (16) and (17), one gets where er f −1 (x) is the inverse function to the error integral er f (x). Assuming that the material initially is free of damage, i.e., when P(0) = 1, one can simplify the latter equation.
Discounting the very tail of the initial crack distribution, one gets max V l 0 ≤ l c . Thus, for well-developed cracks and, in many cases, even for small cracks, the second term in (5) for l 0c dominates the first one. It means that the dependence of l 0c on l c and, therefore, on the material fracture toughness K f can be neglected. Then equation (18) can be approximated by Taking into account that in the case of contact fatigue k 10 is proportional to the maximum contact pressure q max and that it also depends on the friction coefficient λ and the ratio of residual stress q 0 and q max = p H as well as taking into account the relationships μ ln = (where μ and σ are the regular initial mean and standard deviation) one arrives at a simple analytical formula where C 0 depends only on the friction coefficient λ and the ratio of the residual stress q 0 and the maximum Hertzian pressure p H . Finally, assuming that σ μ from (20) one can obtain the formula Also, formulas (20) and (21) can be represented in the form of the Lundberg-Palmgren formula (see [1] and the discussion there).
Formula (21) demonstrates the intuitively obvious fact that the fatigue life N is inversely proportional to the value of the parameter g 0 that characterizes the material crack propagation resistance. Equation (21) exhibits a usual for roller and ball bearings as well as for gears dependence of the fatigue life N on the maximum Hertzian pressure p H . Thus, from the well-known experimental data for bearings the range of n values is 20/3 ≤ n ≤ 9. Keeping in mind that usually σ μ, for these values of n contact fatigue life N is practically inverse proportional to a positive power of the mean crack size, i.e. to μ n 2 −1 . Therefore, fatigue life N is a decreasing function of the initial mean crack (inclusion) size μ. This conclusion is valid for any material survival probability P * and is supported by the experimental data discussed in [1]. In particular, ln N is practically a linear function of ln μ with a negative slope 1 − n 2 which is in excellent agreement with the Timken Company test data [9]. Keeping in mind that n > 2, at early stages of fatigue failure, i.e. when er f −1 (2P * − 1) > 0 for P * > 0.5, one easily determines that fatigue life N is a decreasing function of the initial standard deviation of crack sizes σ. Similarly, at late stages of fatigue failure, i.e. when P * < 0.5, the fatigue life N is an increasing function of the initial standard deviation of crack sizes σ. According to (21), for P * = 0.5 fatigue life N is independent from σ, however, according to (20), for P * = 0.5 fatigue life N is a slowly increasing function of σ. By differentiating p m (N) obtained from (16) with respect to σ, one can conclude that the dispersion of P(N) increases with σ.
The stress intensity factor k 1 decreases as the magnitude of the compressive residual stress q 0 increases and/or the magnitude of the friction coefficient λ decreases. Therefore, in (20) and 21) the value of C 0 is a monotonically decreasing function of residual stress q 0 and friction coefficient λ.
Being applied to bearings and/or gears the described statistical contact fatigue model can be used as a research and/or engineering tool in pitting modeling. In the latter case, some of the model parameters may be assigned certain fixed values based on the scrupulous analysis of steel quality and quality and stability of gear and bearing manufacturing processes.
In case of structural fatigue the Hertzian stress in formulas (20) and (21) should be replaced the dominant stress acting on the part while constant C 0 would be dependent on the ratios of other external stresses acting on the part at hand to the dominant stress in a certain way (see examples of torsional and bending fatigue below).

Examples of torsional and bending fatigue
Suppose that in a beam material the defect distribution is space-wise uniform and follows equation (1). Also, let us assume that the residual stress is zero.
First, let us consider torsional fatigue. Suppose a beam is made of an elastic material with elliptical cross section (a and b are the ellipse semi-axes, b < a) and directed along the y-axis. The beam is under action of torque M y about the y-axis applied to its ends. The side surfaces of the beam are free of stresses. Then it can be shown (see Lurye [10], p. 398) that where G is the material shear elastic modulus, G = E/[2(1 + ν)] (E and ν are Young's modulus and Poisson's ratio of the beam material), and γ is a dimensionless constant. By introducing the principal stresses σ 1 , σ 2 , and σ 3 that satisfy the equation For the case of a > b the maximum principal tensile stress σ 1 = − 2Gγa 2 b a 2 +b 2 is reached at the surface of the beam at points (0, y, ±b) and depending on the sign of M y it acts in one of the directions described by the directional cosines where α is the direction along one of the principal stress axes. For the considered case of elliptic beam, the moments of inertia of the beam elliptic cross section about the x-and y-axes, I x and I z as well as the moment of torsion M y applied to the beam are as follows (see Lurye [10], pp. 395, 399) . Keeping in mind that according to Hasebe and Inohara [11] and Isida [12], the stress intensity factor k 1 for an edge crack of radius l and inclined to the surface of a half-plane at the angle of π/4 (see Then, fatigue life of a beam under torsion follows from substituting the expression for k 10 Now, let us consider bending fatigue of a beam/console made of an elastic material with elliptical cross section (a and b are the ellipse semi-axes) and length L. The beam is directed along the y-axis and it is under the action of a bending force P x directed along the x-axis which is applied to its free end. The side surfaces of the beam are free of stresses. The other end y = 0 of the beam is fixed. Then it can be shown (see Lurye [10]) that where I z is the moment of inertia of the beam cross section about the z-axis. By introducing the principal stresses that satisfy the equation σ 3 − σ y σ 2 − (τ 2 xy + τ 2 zy )σ = 0, one can find that The tensile principal stress σ 1 reaches its maximum 4|P x |L πa 2 b at the surface of the beam at one of the points (±a, 0, 0) (depending on the sign of load P x ) and is acting along the y-axis -the axis of the beam. Based on equations (28) and the solution for the surface crack inclined to the surface of the half-space at angle of π/2 (see Hasebe and Inohara [11] and Isida [12]), one . Therefore, bending fatigue life of a beam follows from substituting the expression for k 10 where function g(μ, σ) is determined by equation (26).
In both cases of torsion and bending, fatigue life is independent of the elastic characteristic of the beam material (see formulas (25), (29), and (26)), and it is dependent on fatigue parameters of the beam material (n and g 0 ), the initial defect distribution (i.e. on μ and σ), the geometry of the beam cross section (a and b), and its length L and the applied loading (P x or M y ).
In a similar fashion the model can be applied to contact fatigue if the stress field is known. A more detailed analysis of contact fatigue is presented below for a two-dimensional case.

Contact problem for an elastic half-plane weakened by straight cracks
A general theory of a stress state in an elastic plane with multiple cracks was proposed in [13].
In this section this theory is extended to the case of an elastic half-plane loaded by contact and residual stresses [1]. A study of lubricant-surface crack interaction, a discussion of the difference between contact fatigue lives of drivers and followers, the surface and subsurface initiated fatigue as well as fatigue of rough surfaces can be found in [1].
The main purpose of the section is to present formulations for the contact and fracture mechanics problems for an elastic half-plane weakened by subsurface cracks. The problems for surface cracks in an elastic lubricated half-plane are formulated and analyzed in [1]. The problems are reduced to systems of integro-differential equations with nonlinear boundary conditions in the form of alternating equations and inequalities. An asymptotic (perturbation) method for the case of small cracks is applied to solution of the problem and some numerical examples for small cracks are presented.
Let us introduce a global coordinate system with the x 0 -axis directed along the half-plane boundary and the y 0 -axis perpendicular to the half-plane boundary and pointed in the direction outside the material. The half-plane occupies the area of y 0 ≤ 0. Let us consider a contact problem for a rigid indenter with the bottom of shape y 0 = f (x 0 ) pressed into the elastic half-plane (see Fig. 2). The elastic half-plane with effective elastic modulus , E and ν are the half-plane Young's modulus and Poisson's ratio) is weakened by N straight cracks. The crack faces are frictionless. Besides the global coordinate system we will introduce local orthogonal coordinate systems for each straight crack of half-length l k in such a way that their origins are located at the crack centers with complex coordinates z 0 k = x 0 k + iy 0 k , k = 1, . . . , N, the x k -axes are directed along the crack faces and the y k -axes are directed perpendicular to them. The cracks are inclined to the positive direction of the x 0 -axis at the angles α k , k = 1, . . . , N. All cracks are considered to be subsurface. The faces of every crack may be in partial or full contact with each other. The indenter is loaded by a normal force P and may be in direct contact with the half-plane or separated from it by a layer of lubricant. The indenter creates a pressure p(x 0 ) and frictional stress τ(x 0 ) distributions. The frictional stress τ(x 0 ) between the indenter and the boundary of the half-plane is determined by the contact pressure p(x 0 ) through a certain relationship. The cases of dry and fluid frictional stress τ(x 0 ) are considered in [1]. At infinity the half-plane is loaded by a tensile or compressive (residual) stress σ ∞ x 0 = q 0 which is directed along the x 0 -axis. In this formulation the problem is considered in [1].
Then the problem is reduced to determining of the cracks behavior. Therefore, in dimensionless variables the equations of the latter problem for an elastic half-plane weakened by cracks and loaded by contact and residual stresses have the following form where the kernels in these equations are described by formulas k , X n = x n e iα n + z 0 n , k, n = 1, . . . , N, where v k (x k ), u k (x k ), and p nk (x k ), k = 1, . . . , N, are the jumps of the normal and tangential crack face displacements and the normal stress applied to crack faces, respectively, a and b are the dimensionless contact boundaries, δ k is the dimensionless crack half-length, δ k = l k /b, δ nk is the Kronecker tensor (δ nk = 0 for n = k, δ nk = 1 for n = k), i is the imaginary unit, i = √ −1.
For simplicity primes at the dimensionless variables are omitted. The characteristic values q andb that are used for scaling are the maximum Hertzian pressure p H and the Hertzian contact half-width a H where R can be taken as the indenter curvature radius at the center of its bottom.
To simplify the problem formulation it is assumed that for small subsurface cracks (i.e. for δ 0 1, δ 0 = max 1≤k≤N δ k ) the pressure p(x 0 ) and frictional stress τ(x 0 ) are known and are close to the ones in a contact of this indenter with an elastic half-plane without cracks. It is worth mentioning that cracks affect the contact boundaries a and b and the pressure distribution p(x 0 ) as well as each other starting with the terms of the order of δ 0 1.
Therefore, for the given shape of the indenter f (x 0 ), pressure p(x 0 ), frictional stress functions τ(x 0 ), residual stress q 0 , crack orientation angles α k and sizes δ k , and the crack positions z 0 k , k = 1, . . . , N, the solution of the problem is represented by crack faces displacement jumps u k (x k ), v k (x k ), and the normal contact stress p nk (x k ) applied to the crack faces (k = 1, . . . , N). After the solution of the problem has been obtained, the dimensionless stress intensity factors k ± 1k and k ± 2k are determined according to formulas

Problem solution
Solution of this problem is associated with formidable difficulties represented by the nonlinearities caused by the presence of the free boundaries of the crack contact intervals and the interaction between different cracks. Under the general conditions solution of this problem can be done only numerically. However, the problem can be effectively solved with the use of just analytical methods in the case when all cracks are small in comparison with the characteristic sizeb of the contact region, i.e., when δ 0 = max 1≤k≤N δ k 1. In this case, it can be shown that the influence of the presence of cracks on the contact pressure is of the order of O(δ 0 ) and with the precision of O(δ 0 ) the crack system in the half-plane is subjected to the action of the contact pressure p 0 (x 0 ) and frictional stress τ 0 (x 0 ) that are obtained in the absence of cracks. The further simplification of the problem is achieved under the assumption that cracks are small in comparison to the distances between them, i.e.
The latter assumption with the precision of O(δ 2 0 ), δ 0 1, provides the conditions for considering each crack as a single crack in an elastic half-plane while the crack faces are loaded by certain stresses related to the contact pressure p 0 (x 0 ), contact frictional stress τ 0 (x 0 ), and the residual stress q 0 . The crucial assumption for simple and effective analytical solution of the considered problem is the assumption that all cracks are subsurface and much smaller in size than their distances to the half-plane surface Essentially, that assumption permits to consider each crack as a single crack in a plane (not a half-plane) with faces loaded by certain stresses related to p 0 (x 0 ), τ 0 (x 0 ), and q 0 .
Let us assume that the frictional stress τ 0 (x 0 ) is determined by the Coulomb law of dry friction which in dimensionless variables can be represented by where λ is the coefficient of friction. In (38) we assume that λ ≥ 0, and, therefore, the frictional stress is directed to the left. It is well known that for small friction coefficients λ the distribution of pressure is very close to the one in a Hertzian frictionless contact. Therefore, the expression for the pressure p 0 (x 0 ) in the absence of cracks with high accuracy can be taken in the form Let us consider the process of solution of the pure fracture mechanics problem described by equations (31)-(33), (35), (38), (39). For small cracks, i.e. for δ 0 1, the kernels from (32) and (33 ) are regular functions of t, x n , and x k and they can be represented by power series in δ k 1 and δ n 1 as follows In (40) where functions v kj , u kj , p nkj have to be determined in the process of solution. Expanding the terms of the equations (31)-(33), and (35) and equating the terms with the same powers of δ 0 we get a system of boundary-value problems for integro-differential equations of the first kind which can be easily solved by classical methods [1,14]. We will limit ourselves to determining only the first two terms of the expansions in (42) in the case of Coulomb's friction law given by (38) and (39). Without getting into the details of the solution process (which can be found in [1]) for the stress intensity factors k ± 1n and k ± 2n we obtain the following analytical formulas [1] k where the kernels are determined according to the formulas and θ(x) is the step function (θ(x) = −1 for x < 0 and θ(x) = 1 for x ≥ 0).

Comparison of analytical asymptotic and numerical solutions for small subsurface cracks
Let us compare the asymptotically (k ± 1a and k ± 2a ) and numerically (k ± 1n and k ± 2n ) obtained solutions of the problem for the case when y 0 = −0.4, δ 0 = 0.1, α = π/2, λ = 0.1, and q 0 = −0.005. The numerical method used for calculating k ± 1n and k ± 2n is described in detail in [1]. Both the numerical and asymptotic solutions are represented in Fig. 3. It follows from Fig. 3 that the asymptotic and numerical solutions are almost identical except for the region where the numerically obtained k + 1 (x 0 ) is close to zero. The difference is mostly caused by the fact that the used asymptotic solution involve only two terms, i.e., the accuracy of these asymptotic solutions is O(δ 2 0 ) for small δ 0 . However, according to the two-term asymptotic solutions the maximum values of k ± 1 differ from the numerical ones by no more than 1.4%. One can expect to get much higher precision if δ 0 < 0.1 and | y 0 | δ 0 .

Stress intensity factors k ± 1n and k ± 2n behavior for subsurface cracks
Some examples of the behavior of the stress intensity factors k ± 1n and k ± 2n for subsurface cracks are presented below.
It is important to keep in mind that for the cases of no friction (λ = 0) and compressive or zero residual stress (q 0 ≤ 0) all subsurface cracks are closed and, therefore, at their tips k ± 1n = 0. Let us consider the case when the residual stress q 0 is different from zero. The residual stress influence on k + 1n results in increase of k + 1n for a tensile residual stress q 0 > 0 or its decrease for a compressive residual stress q 0 < 0 of the material region with tensile stresses. From formulas (43), (44), and Fig. 4 (obtained for y 0 n = −0.2, α n = π/2, and δ n = 0.1) follows that for all x 0 n and for increasing residual stress q 0 (see the curves marked with 3 and 5 that correspond to λ = 0.1, q 0 = 0.04, and λ = 0.2, q 0 = 0.02, respectively) the stress intensity factor k + 1n is a non-decreasing function of q 0 . Moreover, if at some material point k + 1n (q 0 1 ) > 0 for some residual stress q 0 1 , then k + 1n (q 0 2 ) > k + 1n (q 0 1 ) for q 0 2 > q 0 1 (compare curves marked with 1 and 2 with curves marked with 3 and 4 as well as with curves marked with 5 and 6, respectively). Similarly, for all x 0 n when the magnitude of the compressive residual stress (q 0 < 0) increases (see curves marked with 4 and 6 that correspond to λ = 0.1, q 0 = −0.01 and  [5]). Reprinted with permission of Springer. λ = 0.2, q 0 = −0.03, respectively) and at some material point k + 1n (q 0 1 ) > 0 for some residual stress q 0 1 then k + 1n (q 0 2 ) < k + 1n (q 0 1 ) for q 0 2 < q 0 1 . Based on the fact that the normal stress intensity factor k 0+ 1n is positive only in the near surface material layer, we can make a conclusion that this layer increases in size and, starting with a certain value of tensile residual stress q 0 > 0, crack propagation becomes possible at any depth beneath the half-plane surface. For increasing compressive residual stresses, the thickness of the material layer where k + 1n > 0 decreases. The analysis of the results for subsurface cracks following from formulas (43) and (44) shows [1] that the values of the stress intensity factors k ± 1n and k ± 2n are insensitive to even relatively large variations in the behavior of the distributions of the pressure p(x 0 ) and frictional stress τ(x 0 ). In particular, the stress intensity factors k ± 1n and k ± 2n for the cases of dry and fluid (lubricant) friction as well as for the cases of constant pressure and frictional stress are very close to each other as long as the normal force (integral of p(x 0 ) over the contact region) and the friction force (integral of τ over the contact region) applied to the surface of the half-plane are the same [1].
Qualitatively, the behavior of the normal stress intensity factors k ± 1n for different angles of orientation α n is very similar while quantitatively it is very different. An example of that is presented for a horizontal (α n = 0) subsurface crack in Fig. 5 and for a subsurface crack perpendicular to the half-plane boundary (α n = π/2) in  For the same loading and crack parameters the behavior of the shear stress intensity factor k ± 2n is represented in Fig. 7 and 8. It is important to observe that the shear stress intensity factors k ± 2n are insensitive to changes of the friction coefficient λ. Fig. 5 and 6 clearly show that for subsurface cracks with angle α n = π/2 for zero or tensile residual stress q 0 the normal stress intensity factors k ± 1n are significantly higher (by two orders of magnitude) than the ones for α n = 0. For α n = 0 and α n = π/2 the orders of magnitude of the shear stress intensity factors k ± 2n are the same (see Fig. 7 and 8). Moreover, from these graphs it is clear that the normal stress intensity factors k ± 1n are significantly influenced by the friction coefficient λ while the shear stress intensity factors k ± 2n are insensitive to the value of the friction coefficient λ.
Obviously, in the single-term approximation the behavior of k 0− 1n and k 0− 2n is identical to the one of k 0+ 1n and k 0+ 2n , respectively. Generally, the difference between k 0− 1n and k 0+ 1n as well as between k 0− 2n and k 0+ 2n is of the order of magnitude of δ 0 1.

Behavior
The process of lubricant-surface crack interaction is very complex and the details of the problem formulation, the numerical solution approach, and a comprehensive analysis of the results can be found in [1]. Therefore, here we will discuss only the most important features of this phenomenon. It is well known that the presence of lubricant between surfaces in contact is very beneficial as it reduces the contact friction and wear and facilitates better heat transfer from the contact. However, in some cases the lubricant presence may play a detrimental role. In particular, in cases when the elastic solid (half-plane) has a surface crack inclined toward the incoming high contact pressure transmitted through the lubricant. Such a crack may open up and experience high lubricant pressure applied to its faces. This pressure creates the normal stress intensity factor k − 1n far exceeding the value of the stress intensity factors k ± 1n for comparable subsurface cracks while the shear stress intensity factors for surface k − 2n and and subsurface k ± 2n cracks remain comparable in value. That becomes obvious from the comparison of the graphs from Fig. 4-8 with the graphs from Fig. 10 and 9.
In cases when a surface crack is inclined away from the incoming high lubricant pressure the crack does not open up toward the incoming lubricant with high pressure and it behaves similar to a corresponding subsurface crack, i.e. its normal stress intensity factor k − 1n in its value is similar to the one for a corresponding subsurface crack. It is customary to see the normal stress intensity factor for surface cracks which open up toward the incoming high lubricant pressure to exceed the one for comparable subsurface cracks by two orders of magnitude. In such cases the compressive residual stress has very little influence on the crack behavior due to domination of the lubricant pressure.
The possibility of high normal stress intensity factors for surface cracks leads to serious consequences. In particular, it explains why fatigue life of drivers is usually significantly higher than the one for followers [1]. Also, it explains why under normal circumstances fatigue failure is of subsurface origin [1].

Stress intensity factors and directions of fatigue crack propagation
The process of fatigue failure is usually subdivided into three major stages: the nucleation period, the period of slow pre-critical fatigue crack growth, and the short period of fast unstable crack growth ending in material loosing its integrity. The durations of the first two stages of fatigue failure depend on a number of parameters such as material properties, specific environment, stress state, temperature, etc. Usually, the nucleation period is short [1]. We are interested in the main part of the process of fatigue failure which is due to slow pre-critical crack growth. In these cases max (k + 1 , k − 1 ) < K f , where K f is the material fracture toughness.
During the pre-critical fatigue crack growth cracks remain small. Therefore, they can be modeled by small straight cuts in the material. For such small subsurface cracks it is sufficient to use the one-term asymptotic approximations k ± 1 = k 1 and k ± 2 = k 2 for the stress intensity factors from (43) and (44) which in dimensional variables take the form where i is the imaginary unit (i 2 = −1), θ(x) is a step function: θ(x) = 0, x ≤ 0 and θ(x) = 1, x > 0. It is important to mention that according to (45) for subsurface cracks the quantities of k 10 = k 1 l −1/2 and k 20 = k 2 l −1/2 are functions of x and y and are independent from l.
Numerous experimental studies have established the fact that at relatively low cyclic loads materials undergo the process of pre-critical failure while the rate of crack growth dl/dN (N is the number of loading cycles) in the predetermined direction is dependent on k ± 1 and K f . A number of such equations of pre-critical crack growth and their analysis are presented in [6]. However, what remains to be determined is the direction of fatigue crack growth.
Assuming that fatigue cracks growth is driven by the maximum principal tensile stress (see the section on Three-Dimensional Model of Contact and Structural Fatigue) we immediately obtain the equation which determines the orientation angles α ± of a fatigue crack growth at the crack tips. Due to the fact that a fatigue crack remains small during its pre-critical growth (i.e. practically during its entire life span) and being originally modeled by a straight cut with half-length l at the point with coordinates (x, y) the crack remains straight, i.e. the crack direction is characterized by one angle α = α + = α − . This angle is practically independent from crack half-length l because k ± 2 = k ± 20 √ l, where k ± 20 is almost independent from l for small l (see (45)). The dependence of k ± 2 on the number of loading cycles N comes only through the dependence of the crack half-length l on N. Therefore, the crack angle α is just a function of the crack location (x, y).
Along these directions k 1 reaches its extremum values. The actual direction of crack propagation α is determined by one of these two angles α 1 and α 2 for which the value of the normal stress intensity factor k 1 (N, x, y, l, α) is greater.
A more detailed analysis of the directions of fatigue crack propagation can be found in [1].

Two-dimensional contact fatigue model
In a two-dimensional case compared to a three-dimensional case a more accurate description of the contact fatigue process can be obtained due to the fact that in two dimensions it is relatively easy to get very accurate formulas for the stress intensity factors at crack tips [1,5]. The rest of the fatigue modeling can be done the same way as in the three-dimensional case with few simple changes. In particular, in a two-dimensional case of contact fatigue only subsurface originated fatigue is considered and cracks are modeled by straight cuts with half-length l. That gives the opportunity to use equations (45) Figure 10. Distributions of the stress intensity factors k − 1 (curve 1) and k − 2 (curve 2) in case of a surface crack: α = π/6, λ = 0.1, and q 0 = −0.5 (after Kudish and Burris [16]). Reprinted with permission from Kluwer Academic Publishers. factors k 1 , k 2 , and crack angle orientations α. The rest of contact fatigue modeling follows exactly the derivation presented in the section on Three-Dimensional Model of Contact and Structural Fatigue. Therefore, the fatigue life N with survival probability P * of a contact subjected to cyclic loading can be expressed in the form [1,17] where er f −1 (x) is the inverse function to the error integral er f (x) [8].
First, let us consider the model behavior in some simple cases. If f (0, x, y, z, l 0 ) is a uniform crack distribution over the material volume V (except for a thin surface layer where f = 0). Then based on (11) it can be shown that p(N, x, y, z) reaches its minimum at the points where k 10 and the principal tensile stress reach their maximum values. This leads to the conclusion that the material local failure probability (1 − p) reaches its maximum at the points with maximal tensile stress. Therefore, for a uniform initial crack distribution f (0, x, y, z, l 0 ) the survival probability P(N) from (16) is determined by the material local survival probability at the points at which the maximal tensile stress is attained.
However, the latter conclusion is not necessarily correct if the initial crack distribution f (0, x, y, z, l 0 ) is not uniform over the material volume. Suppose, k 10 (x, y, z) is maximal at the point (x m , y m , z m ) and at the initial time moment N = 0 at some point (x * , y * , z * ) there exist cracks larger than the ones at the point (x m , y m , z m ), namely, Then after a certain number of loading cycles N > 0 the material damage at point (x * , y * , z * ) may be greater than at point (x m , y m , z m ), where l 0c reaches its maximum value. Therefore, fatigue failure may occur at the point (x * , y * , z * ) instead of the point (x m , y m , z m ), and the material weakest point is not necessarily is the material most stressed point.
If μ ln and σ ln depend on the coordinates of the material point (x, y, z), then there may be a series of points where in formula (16) for the given number of loading cycles N the minimum over the material volume V is reached. The coordinates of such points may change with N. This situation represents different potentially competing fatigue mechanisms such as pitting, flaking, etc. The occurrence of fatigue damage at different points in the material depends on the initial defect distribution, applied stresses, residual stress, etc.
In the above model of contact fatigue the stressed volume V plays no explicit role. However, implicitly it does. In fact, the initial crack distribution f (0, x, y, z, l 0 ) depends on the material volume. In general, in a larger volume of material, there is a greater chance to find inclusions/cracks of greater size than in a smaller one. These larger inclusions represent a potential source of pitting and may cause a decrease in the material fatigue life of a larger material volume.
Assuming that μ ln and σ ln are constants, and assuming that the material failure occurs at the point ( where parameter n can be compared with constant c/e in the Lundberg-Palmgren formula [1]. The major difference between the Lundberg-Palmgren formula and formula (49) derived from this model of contact fatigue is the fact that in (49) constant C * depends on material defect parameters μ, σ, coefficient of friction λ, residual stresses q 0 , and probability of survival P * in a certain way while in the Lundberg-Palmgren formula the constant C * depends only on the depth z 0 of the maximum orthogonal stress, stressed volume V, and probability of survival P * .
Let us analyze formula (21). It demonstrates the intuitively obvious fact that the fatigue life N is inverse proportional to the value of the parameter g 0 that characterizes the material crack propagation resistance. So, for materials with lower crack propagation rate, the fatigue life is higher and vice versa. Equation (21) exhibits a usual for gears and roller and ball bearings dependence of fatigue life N on the maximum Hertzian pressure p H . Thus, from the well-known experimental data for bearings, the range of n values is 20/3 ≤ n ≤ 9. Keeping in mind that usually σ μ, for these values of n contact fatigue life N is practically inverse proportional to a positive power of the crack initial mean size, i.e., to μ n/2−1 . Therefore, fatigue life N is a decreasing function of the initial mean crack (inclusion) size μ and ln N = −(n/2 − 1) ln μ + constant. This conclusion is valid for any value of the material survival probability P * and is supported by the experimental data obtained at the Timken Company by Stover and Kolarik II [18] and represented in Fig. 11 in a log-log scale. If P * > 0.5 then er f −1 (2P * − 1) > 0 and (keeping in mind that n > 2) fatigue life N is a decreasing function of the initial standard deviation of crack sizes σ. Similarly, if P * < 0.5, then fatigue life N is an increasing function of the initial standard deviation of crack sizes σ. According to (20), for P * = 0.5 fatigue life N is independent from σ, and, according to (21), for P * = 0.5 fatigue life N is a slowly increasing function of σ. By differentiating p m (N) obtained from (16) with respect to σ, we can conclude that the dispersion of P(N) increases with σ.
From (45) (also see Kudish [1]) follows that the stress intensity factor k 1 decreases as the magnitude of the compressive residual stress q 0 increases and/or the magnitude of the friction coefficient λ decreases. Therefore, it follows from formulas (45) that C 0 is a monotonically decreasing function of the residual stress q 0 and friction coefficient λ. Numerical simulations of fatigue life show that the value of C 0 is very sensitive to the details of the residual stress distribution q 0 versus depth.
Let us choose a basic set of model parameters typical for bearing testing: maximum Hertzian pressure p H = 2 GPa, contact region half-width in the direction of motion a H = 0.249 mm, friction coefficient λ = 0.002, residual stress varying from q 0 = −237.9 MPa on the surface to q 0 = 0.035 MPa at the depth of 400 μm below it, fracture toughness K f varying between 15 and 95 MPa · m 1/2 , g 0 = 8.863 MPa −n · m 1−n/2 · cycle −1 , n = 6.67, mean of crack initial half-lengths μ = 49.41 μm (μ ln = 3.888 + ln(μm)), crack initial standard deviation σ = 7.61 μm (σ ln = 0.1531). Numerical results show that the fatigue life is practically independent from the material fracture toughness K f , which supports the assumption used for the derivation of formulas Pitting probability 1 − P(N) calculated for the basic set of parameters (solid curve) with μ = 49.41 μm, σ = 7.61 μm (μ ln = 3.888 + ln(μm), σ ln = 0.1531), for the same set of parameters and the increased initial value of crack mean half-lengths (dash-dotted curve) μ = 74.12 μm (μ ln = 4.300 + ln(μm), σ ln = 0.1024), and for the same set of parameters and the increased initial value of crack standard deviation (dotted curve) σ = 11.423 μm (μ ln = 3.874 + ln(μm), σ ln = 0.2282) (after Kudish [17]). Reprinted with permission from the STLE. Pitting probability 1 − P(N) calculated for the basic set of parameters (solid curve) and for the same set of parameters and changed profile of residual stress q 0 (dashed curve) in such a way that at points where q 0 is compressive its magnitude is unchanged and at points where q 0 is tensile its magnitude is doubled (after Kudish [17]). Reprinted with permission from the STLE. parameters, just one parameter from the basic set of parameters will be varied at a time and graphs of the pitting probability 1 − P(N) for these sets of parameters (basic and modified) will be compared. Figure 12 shows that as the initial values of the mean μ of crack half-lengths and crack standard deviation σ increase contact fatigue life N decreases. Similarly, contact fatigue life decreases as the magnitude of the tensile residual stress and/or friction coefficient increase (see Fig. 4 and 13). The results show that the fatigue life does not change when the magnitude of the compressive residual stress is increased/decreased by 20% of its base value while the tensile portion of the residual stress distribution remains the same. Obviously, that is in agreement with the fact that tensile stresses control fatigue. Moreover, the fatigue damage occurs in the region with the resultant tensile stresses close to the boundary between tensile and compressive residual stresses. However, when the compressive residual stress becomes small enough the acting frictional stress may supersede it and create new regions with tensile stresses that potentially may cause acceleration of fatigue failure.  Table 1. Relationship between the tapered bearing fatigue life N 15.9 and the initial inclusion size mean and standard deviation (after Kudish [17]). Reprinted with permission from the STLE.
Let us consider an example of the further validation of the new contact fatigue model for tapered roller bearings based on a series of approximate calculations of fatigue life. The main simplifying assumption made is that bearing fatigue life can be closely approximated by taking into account only the most loaded contact. The following parameters have been used for calculations: p H = 2.12 GPa, a H = 0.265 mm, λ = 0.002, g 0 = 6.009 MPa −n · m 1−n/2 · cycle −1 , n = 6.67, the residual stress varied from q 0 = −237.9 MPa on the surface to q 0 = 0.035 MPa at the depth of 400 μm below the surface, fracture toughness K f varied between 15 and 95MPa · m 1/2 . The crack/inclusion initial mean half-length μ varied between 49.41 and 244.25 μm (μ ln = 3.888 − 5.498 + ln(μm)), the crack initial standard deviation varied between σ = 7.61 and 37.61 μm (σ ln = 0.1531). The results for fatigue life N 15.9 (for P(N 15.9 ) = P * = 0.159) calculations are given in the Table 1 and practically coincide with the experimental data obtained by The Timken Company and presented in Fig. 19 by Stover, Kolarik II, and Keener [? ] (in the present text this graph is given as Fig. 11). One must keep in mind that there are certain differences in the numerically obtained data and the data presented in the above mentioned Fig. 11 due to the fact that in Fig. 11 fatigue life is given as a function of the cumulative inclusion length (sum of all inclusion lengths over a cubic inch of steel) while in the model fatigue life is calculated as a function of the mean inclusion length.
It is also interesting to point out that based on the results following from the new model, bearing fatigue life can be significantly improved for steels with the same cumulative inclusion length but smaller mean half-length μ (see Fig. 12). In other words, fatigue life of a bearing made from steel with large number of small inclusions is higher than of the one made of steel with small number of larger inclusions given that the cumulative inclusion length is the same in both cases. Moreover, bearing and gear fatigue lives with small percentage of failures (survival probability P > 0.5) for steels with the same cumulative inclusion length can also be improved several times if the width of the initial inclusion distribution is reduced, i.e., when the standard deviation σ of the initial inclusion distribution is made smaller (see Figure 12). Figures 13 and 14 show that the elevated values of the tensile residual stress are much more detrimental to fatigue life than greater values of the friction coefficient.
Finally, the described model is flexible enough to allow for replacement of the density of the initial crack distribution (see (1)) by a different function and of Paris's equation for fatigue crack propagation (see (4)) by another equation. Such modifications would lead to results on fatigue life varying from the presented above. However, the methodology, i.e., the way the formulas for fatigue life are obtain and the most important conclusions will remain the same.
This methodology has been extended on the cases of non-steady cyclic loading as well as on the case of contact fatigue of rough surfaces [1]. Also, this kind of modeling approach has been applied to the analysis of wear and contact fatigue in cases of lubricant contaminated by rigid abrasive particles and contact surfaces charged with abrasive particles [19] as well as to calculation of bearing wear and contact fatigue life [20].
In conclusion we can state that the presented statistical contact and structural fatigue models take into account the most important parameters of the contact fatigue phenomenon (such as normal and frictional contact and residual stresses, initial statistical defect distribution, orientation of fatigue crack propagation, material fatigue resistance, etc.). The models allows for examination of the effect of variables such as steel cleanliness, applied stresses, residual stress, etc. on contact fatigue life as single or composite entities. Some analytical results illustrating these models and their validation by the experimentally obtained fatigue life data for tapered bearings are presented.

Closure
The chapter presents a detailed analysis of a number of plane crack mechanics problems for loaded elastic half-plane weakened by a system of cracks. Surface and subsurface cracks are considered. All cracks are considered to be straight cuts. The problems are analyzed by the regular asymptotic method and numerical methods. Solutions of the problems include the stress intensity factors. The regular asymptotic method is applied under the assumption that cracks are far from each other and from the boundary of the elastic solid. It is shown that the results obtained for subsurface cracks based on asymptotic expansions and numerical solutions are in very good agreement. The influence of the normal and tangential contact stresses applied to the boundary of a half-plane as well as the residual stress on the stress intensity factors for subsurface cracks is analyzed. It is determined that the frictional and residual stresses provide a significant if not the predominant contribution to the problem solution. Based on the numerical solution of the problem for surface cracks in the presence of lubricant the physical nature of the "wedge effect" (when lubricant under a sufficiently large pressure penetrates a surface crack and ruptures it) is considered. Solution of this problem also provides the basis for the understanding of fatigue crack origination site (surface versus subsurface) and the difference of fatigue lives of drivers and followers. New twoand three-dimensional models of contact and structural fatigue are developed. These models take into account the initial crack distribution, fatigue properties of the solids, and growth of fatigue cracks under the properly determined combination of normal and tangential contact and residual stresses. The formula for fatigue life based on these models can be reduced to a simple formula which takes into account most of the significant parameters affecting contact and structural fatigue. The properties of these contact fatigue models are analyzed and the results based on them are compared to the experimentally obtained results on contact fatigue for tapered bearings.