COMPUTATIONAL FATIGUE ANALYSIS OF CONTACTING MECHANICAL ELEMENTS

Original scientific paper The paper investigates the influence of different parameters on stress distribution in the contact area of mechanical elements. The study is aimed to determine particular positions in the contacting bodies, where small crack should be initiated due to large stresses resulting from different contact conditions. A two-dimensional equivalent contact model together with Hertzian contact theory is used for all computational analyses utilising the finite element method. The distribution of the maximum equivalent von Mises stress under the contacting area is first determined for pure Hertzian contact conditions, where the maximum stresses appear at a certain depth under the contacting surface. Then a parametric study to determine the influence of loading, contact surface curvature, contact friction, residual stresses and elasto-hydro-dynamic lubrication on the stress distribution is performed. They tend to increase the stresses and move their maximum to or very close to the contacting surface, which can lead to sub-surface or surface crack initiation. In this study the surface crack initiation and propagation on the gear teeth flank is presented. The virtual crack extension (VCE) method, implemented in the finite element method, is used for simulating the fatigue crack propagation from the initial crack up to the formation of the surface pit. Comparison between the results of the numerical simulation and the results of experimental testing shows good correlation.


Introduction
Contact mechanics is concerned with stresses and deformations that arise when the surfaces of two solid bodies are brought into contact [1].The contact can be either conforming or non-conforming.Conforming contact is the contact where the surfaces of the two bodies fit each other exactly or very closely without deformation.In the contrast, the contact of bodies with dissimilar profiles is characterised as non-conforming.The contact area of the non-conforming contact is generally small comparing with the dimensions of the bodies themselves.The stresses are highly concentrated in the region close to the contact zone and are not greatly influenced by the shape of the bodies at a distance from the contact area.Common engineering applications of non-conforming contacts are wheels, bearings, traction drives, gears and cams [2,3].
The maximum stresses in contact area resulting from the contact conditions should always be lower than the local material yield stress to prevent any irrecoverable damage of contacting bodies.The material yield stress is usually obtained from a simple one-dimensional tension or shear experimental testing.The load at which plastic yielding begins in realistic structures is related to the uniaxial yield stress through the appropriate yield criterion.For most metallic materials the von Mises criterion is the best choice and is postulated as [4] ( ) ( ) ( ) where σ 1 , σ 2 and σ 3 are the principal stresses, and σ y and τ y are the yield stress of the material in simple tension (or compression) and shear, respectively.The contact conditions between two contacting bodies are complex and many parameters may influence the equivalent stress distribution σ ME in the contact area [5].
This paper reports on the results of computational analyses aimed to determine the influence of loading, contact surface curvature, contact friction, residual stresses and lubrication conditions on distribution of the equivalent stress in the contact area of mechanical elements.The influences of some typical contact parameters, i.e. the maximum contact pressure, the contact surface curvature, coefficient of friction, tensile and compressive residual stresses, mean velocity of contacting surfaces, lubricant viscosity, on the contact stress distribution are clearly illustrated.Such conditions commonly lead to the development of sub-surface or surface initial cracks.
The paper describes a computational simulation of the crack growth leading to pitting which starts from the initial surface crack.The proposed model is based on a 2dimensional plain strain crack propagation analysis, where the required material properties are obtained from common fatigue tests.Although the pitting process involves small amounts of plastic flow below the contact, the crack growth model is used without consideration of friction and plastic flow occurring on the developing crack faces.The prime aim of the study is determination of the fatigue crack growth under various influential contact loading parameters for the case of surface initiated fatigue crack, and the resulting sizes of the estimated pit shapes.As a spring-off, the evaluated functional relationship between the stress intensity factor and the crack length can be used for computational service life predictions of gears in regard to pitting, if coupled with particular short fatigue crack growth theories [6 ÷ 8].

Equivalent contact model
For computational simulations of non-conforming contact of mechanical elements it is advantageous to use a substitute model of two contacting cylinders instead of simulating the actual contact, Fig. 1.The substitute cylinders have the same radii, as are the curvature radii of the treated mechanical elements at the point of the actual contact.The substitute contact model of two cylinders can be transformed into equivalent contact model by utilising the Hertzian contact theory [4].The equivalent contact model consists of a simple cylinder of equivalent radius and material properties along with applied standard Hertzian contact conditions.

Real contact
Hertz idealised contact Following the Hertz contact theory, the distribution of normal pressure in the frictionless contact between two contacting bodies can be analytically determined with [4] , π in which F N is the normal force per unit length of equivalent cylinders and b is the half-width of the contact area, which is given as , π where R * is the equivalent radius and E * is the equivalent Young's modulus defined as .
Various contacting mechanical elements can then be simply simulated by analysing the equivalent cylinder of radius R * and Young's modulus E * (see Fig. 2) with applied contact conditions along the contact area 2b.Such an equivalent problem can be solved with any generalpurpose linear elastic finite element code.This avoids the need to use special algorithms for simulation of the actual contact state, which are not so generally available.The stress distribution under the contact area is fully consistent with the Hertzian analytical results and those obtained using the special contact elements and appropriate solution algorithms [6,7].This leads to the conclusion that using the equivalent contact model shown in Fig. 2 provides adequate results and can therefore be used for determination of the influence of different contact parameters on the contact stress field.3), ( 4) and (5) this results in the half-width of the contact area b = 0,274 mm, the equivalent radius R * = 10 mm and Young's modulus E * = 2,264×10 5 MPa.The used model discretisation is illustrated in Fig. 2, where the finite element mesh has been refined in the contact area to account for large strain gradients appearing in this area.In all computations the plane strain case has been assumed.

Influence of loading and contact surface curvature
Using the simplified computational model presented in Section 2 together with set reference data, the maximum von Mises equivalent stress max ME σ and its position, i.e. depth H under the contact surface, have been computed for different maximum contact pressures p o and equivalent curvature radii R * that are commonly encountered in engineering applications.The results are summarised in Tab. 1.
Computational analyses confirmed the theoretical and experimental findings that the maximum stresses appear at depth H≈ 0,7⋅b in this type of contact.Analysing the results in Table 1, it follows that the maximum von Mises equivalent stress increases with the increase of the maximum contact pressure, while its position under the contact surface depends both on the maximum contact pressure and the equivalent curvature radius.If the maximum contact pressure or the equivalent curvature radius is increased the maximum von Mises stress appears deeper under the contact surface.

Influence of contact friction
In addition to normal contact pressure the contact surfaces of mechanical elements usually transmit a certain tangential traction forces due to friction.The influence of these frictional forces on contact stress can be simulated in its simplest form by applying the Coulomb friction law, where the distribution of frictional (tangential) contact loading q(x) can be easily determined from the normal loading p(x) and coefficient of friction µ with .Simulations have shown that for small values of friction (µ < 0,05) the influence of tangential loading on the equivalent stress σ ME in the contact area is very small and can be neglected.With increase of the coefficient of friction the stress close to the surface also increases.At µ ≈ 0,25 the maximum equivalent stress max ME σ reaches the contact surface and rapidly increases with further increase of µ.

Influence of residual stresses
If the mechanical component is loaded so that it sustains some plastic deformation and is then unloaded, certain residual stresses will remain in the material.For the plane strain conditions the residual stress components r xz τ and r zy τ can be neglected, with remaining stress components becoming independent of x and z axes (see Fig. 1).The equilibrium of residual stresses for a traction free surface then reduces to [4] .0 ), ( , ) (   (9).For ease of application it was assumed that the residual stresses r x σ and r z σ are equal and constant through the surface layer considered.The von Mises equivalent stresses σ EM were then evaluated by Eq. ( 1).The results of these computations for different values of residual stresses are shown in Fig. 4. From Fig. 4 it can be concluded that tension residual stresses are not desirable since they amount to higher equivalent stress σ ME , while compressive residual stresses help to reduce the total equivalent stress σ ME .

Influence of the elasto-hydro-dynamic lubrication
Under the elasto-hydro-dynamic (EHD) lubrication conditions the contacting bodies are separated in a contact zone with a fluid-film, which transmits the contact stresses by its pressure.The transmitting pressure develops in the entry region by hydro-dynamic action and is accompanied by a very large increase in viscosity.The film thickness at the outlet region, i.e. at the end of converging zone, is limited by the necessity of maintaining a finite pressure.This condition virtually determines the film thickness in relation to the mean surface speed u, contact surface radius R and viscous properties of the lubricant [9, 10], see Fig. 5.

Figure 5Conditions in the EHD-lubricated contact
The highly viscous fluid passes through the parallel zone until the pressure and viscosity collapse at the exit, which dictates thinning of the film.The inlet and outlet regions are effectively independent; they meet at the end of the parallel zone with a discontinuity in a contact strain field, which is associated with a sharp pressure peak, i.e. pressure spike.The presence of a pressure spike in the outlet region of the EHD-lubricated contact produces large shear stresses that are localised very close to the surface, and can produce initial cracks on the surface.
The parameters generally used in EHD-lubrication contact analysis are defined as [11 ÷ 13] • Dimensionless load parameter: .
• Dimensionless speed parameter: • Dimensionless material parameter: , where F N is the load per unit contact width, η o is the dynamic viscosity at the atmospheric pressure, u = (u 1 + u 2 )/2 is the mean surface velocity with u 1 and u 2 being respective surface velocities, ρ is the lubricant density, n is the kinematic viscosity and α is the pressure-viscosity coefficient.
Hamrocket al. proposed the following equations for determination of the dimensionless pressure spike amplitude Y and location X [14], see Fig. 5: In all simulations treated in this section it was assumed that the dimensionless load parameter W and the dimensionless material parameter G are constant while the dimensionless speed parameter U varies in terms of different mean surface velocity u and kinematic viscosityn.The lubricant density was equal to ρ = 9×10 −7   kg/m 3 .The distribution of von Mises equivalent stress σ ME under the contact surface has been determined with the computational analysis of the finite element model shown in Fig. 4. The frictionless Hertzian loading has been applied for the maximum contact pressure p o = 1550 MPa with the addition of extra pressure spike due to EHD-lubrication, which was evaluated according to Fig. 5 and Eqs. ( 13) and (14).The influence of the dimensionless speed parameter U on the contact stresses was studied for two cases (see Tab. 2): • Case 1: variable mean velocity and constant kinematic viscosity, • Case 2: constant mean velocity and variable kinematic viscosity.For small values of the mean velocity (u < 10 m/s) its influence on the equivalent stress is very small and can be neglected.With increase of the mean velocity the maximum equivalent stress also increases and its maximum max ME σ moves closer to the contact surface.At u > 30 m/s the maximum equivalent stress max ME σ reaches the contact surface (Fig. 6).The same applies for the influence of the kinematic viscosity.For values of n < 200 mm 2 /s the influence on the stress field can be neglected, while for n > 400 mm 2 /s the maximum equivalent stress moves very close to the surface (Fig. 7).Such conditions commonly lead to the development of surface or subsurface crack initiation.

Initiation and propagation of surface crack
The paper presents a computational model for the simulation of surface initiated fatigue crack growth on gear teeth flanks.A simple contact model is used for simulating fatigue crack growth under conditions of rolling and sliding contact as described in previous sections.The contact model is subjected to moving normal (normal contact pressure) and tangential (frictional forces) contact forces, which also take into account the influence of EHD-lubrication conditions.
The surface-breaking initial cracks are always angled in the opposite direction of the acting frictional forces (see Fig. 8).The finite element mesh shown in Fig. 9 has been used in subsequent analyses.For the configuration of the initial crack on the surface, located at point B, it was assumed that the initial length of the crack is equal to a o =15 µm, with the initial inclination angle towards the contact surface equal to α o = 22°.This configuration follows the metallographic investigations of initial cracks appearing on gears [15].
For the purpose of the crack growth simulation, the virtual crack extension (VCE) method in the framework of the finite element method (FEM) has been applied [16].The VCE method falls into category of energy based methods for crack propagation modelling, where it is assumed that the crack will always propagate in the direction that maximises the release of the strain energy, resulting in the minimum possible total potential energy of the system.Since VCE method does not allow for separation of mode I and mode II stress intensity factors in mixed mode situations, only the combined stress intensity factor K can be determined from the numerically evaluated fracture energy release rate G.The crack is incrementally extended always in the direction of the current maximum potential energy release rate, which necessitates local remeshing around the new crack tip at each increment.The incremental procedure is repeated until full fracture occurs or until the stress intensity factor reaches the critical value K Ic , when full fracture is imminent.
The FE analysis program BERSAFE has been used for computational estimation of the stress intensity factor K and subsequent incremental crack growth simulation.In numerical computations, the crack increment was of size ∆a = 1,5 µm.The stress intensity factor K was estimated in each crack increment for 30 different virtual crack tip extensions (see Fig. 10).Five different loading configurations have been considered in each computation for the purpose of simulating the effect of the moving contact of the gear flanks (see Fig. 11).For each crack increment, the crack was actually extended in the direction of the recorded K max from all calculated load cases.

Figure 12
Stress intensity factor for initial surface micro crack propagation Fig. 12 shows the relationship between stress intensity factor K and the crack length a, and also the shape and magnitude of the surface pit.It can be seen that the computed stress intensity factor K is very small at the beginning, but later increases as the crack propagates towards the contact surface.Numerical simulations have shown that at the moment when the crack reaches the vicinity of the contact surface, the stress intensity factor is extremely high.At that moment it can be expected that the material surface layer breaks away and the pit occurs on the surface.Because of the very small dimensions of surface pits, they can be termed "micro-pitting".Based on experimental testing of the spur gear pair [15,17] with the same operating conditions and loading parameters as used in the numerical computations, the comparison between numerically and experimentally determined pits on the gear flanks was compared.It can be seen that the shape and magnitude of surface pits are in a good agreement [see Fig. 13].

Conclusion
The paper investigates the influences of different parameters on the distribution of stresses in the contact area of mechanical elements.The study is aimed to determine positions of the maximum stresses in the contacting bodies, where the likelihood of irrecoverable damage is the highest.A two-dimensional equivalent Hertzian contact model is used for all computational analyses utilising the finite element method.With aid of computational analyses the parametric study is performed to determine the influence of contact loading, surface curvature, contact friction, residual stresses and elastohydro-dynamic lubrication on the contact stress distribution.From this study it follows that high contact pressures, small contact surface curvatures, high friction, tensile residual stresses, high mean velocities of contacting surfaces and high kinematic viscosity have a negative influence on the stress distribution.They tend to increase the stresses and move their maximum to or very close to the contacting surface, which usually results in development of damage in the form of plastic deformation due to dislocation motion, which can lead to subsurface or surface initial crack development and ultimately to surface pitting.In this paper the surface initial crack and its propagation is investigated in detail.The virtual crack extension (VCE) method in the framework of FEA is used

Figure
Figure 1Transformation from (a) the actual contact to (b) substitute model of two cylinders

Figure 2
Figure 2 Discretised equivalent model for stress field determination in the contact area For all computational simulations treated in this paper a contact of two cylinders with radii R 1 = R 2 = 20 mm, the Young's modulus E 1 = E 2 = 2,06×10 5 MPa and the Poisson's ratio m 1 = m 2 = 0.3 is considered.According to Eqs. (3), (4) and(5) this results in the half-width of the contact area b = 0,274 mm, the equivalent radius R * = 10 mm and Young's modulus E * = 2,264×105 MPa.The used model discretisation is illustrated in Fig.2, where the finite element mesh has been refined in the contact area to account for large strain gradients appearing in this area.In all computations the plane strain case has been assumed.

Figure 3
Figure 3 The influence of the coefficient of friction µ on the equivalent stress σME in the contact area The results of simulations for different values of µ are shown in Fig. 3, where the resulting distribution of the von Mises equivalent stress σ ME along the y-axis of the contact area (depth under the contact surface) is given.In all simulations the maximum contact pressure p o =1550 condition with consideration of residual stresses can be easily evaluated with Eq. (1).

Figure 4
Figure 4 Influence of the residual stresses r x σ and r z σ on the von Mises equivalent stress σME in the contact area To determine the influence of residual stresses on the stress field of the equivalent contact model, the stresses have been computed first by the FEM for the frictionless

Figure 6 Figure 7
Figure 6 Distribution of the von Mises equivalent stress σME in the contact area for kinematic viscosity n= 220 mm 2 /s

Figure 8
Figure 8 Orientation of the initial surface-breaking crack and contact loading

Figure 9
Figure 9 FE discretisation (a) and configuration of the initial crack (b)

Figure 13
Figure 13 Numerically (above) and experimentally determined pit shapes for two-dimensional simulation of fatigue crack propagation from the initial surface crack up to the formation of the surface pit.Comparison of numerically predicted and experimentally recorded pit shapes shows that they are in a very good agreement.

Table 1
The maximum von Mises equivalent stress max