Coupled Numerical Method for Rolling Contact Fatigue Analysis

1560-1567 Abstract: In this paper, our aim was to create a coupled numerical simulation method for rolling contact fatigue (RCF) analysis which takes into account the effects of results of the heat treatment process of the manufacturing of specimen/machine parts, such as hardness, mechanical properties, residual stresses in the carburized layers of the body. Thus, the finite element method (FEM) simulation process calculates the material properties and stress close to the realistic state and using them as an input data, the RCF analysis is carried out by this coupled method. The Dang Van, Matake and Findley multiaxial fatigue criteria were used for failure analysis. In parallel with the development of the method, validation experiments were carried out to decide which criterion is the most appropriate for our purpose and to determine the coupled method efficiency.


INTRODUCTION
The fatigue due to rolling contact is the characteristic damage that determines the life of contact surface pairs, such as bearings or gears, for which the edges of the contact zone are the most critical regions; furthermore damage can appear shortly below the surfaces. Surfaces are mainly subjected to high-cycle or ultrahigh-cycle fatigue which can best be described by the stress-life approach [1].
Multidisciplinary approach is required with application of the most recently developed numerical and testing methods to study the crucial issues of the conditions of the parts and estimate the lifetime [2]. Although the numerical tools and experimental techniques developed a lot in order to be able to characterize the rolling contact events [3], there are still some open issues and gaps that need to be addressed in further research.
Although the answer to questions about the lifetime and condition of the parts is crucial, while numerical and experimental methods are highly evolving and modern testing methods are becoming more and more capable of characterizing such tribological events [4]- [6], there are still problems that require further research. To solve them, it is also necessary to know the mechanical material properties close to the surface, the determination or prediction of which is also a challenge due to the surface treatments. Similarly, knowledge of the residual stress created during the heat treatment is essential as the residual stress can directly modify the fatigue process and thus the lifespan of the examined part [7].
Our goal was to make a coupled numerical simulation method which is not only suitable for determining the real characteristics of the near-surface layer properties and the residual stress, but also for using these results as an initial value to calculate the failure during high-cycle rolling contact. As for model material 16MnCr5 case hardening steel was selected, since it is widely used steel for engineering application. Two-disc rolling fatigue tests for validating the numerical model were carried out with surface treated rollers made from the selected steel.

VALIDATION EXPERIMENTS 2.1 Roller Specimen
The geometry and other properties of the roller specimens was defined based on [8] article. The geometry can be seen in Fig.1. The roller specimens were carburized. The chemical composition of the model material is shown in the Tab.1. The BaCO3 charcoal activator was used as agent for carbonization, the specimens were quenched in oil after carbonization and stress relief annealing was performed. (Fig. 2).

Fatigue Test
Five pairs of rollers were tested with different contact force with the support of Montanuniversität Leoben, Lehrstuhl für Allgemeinen Maschinenbau. One of the tests had to be stopped after one million of rotations due to failure of the shaft. Other tests were running till the acoustic trigger sign of surface damage The clamping force between the discs was 10000 N, 12500 N, 15000 N with p = 2.25 GPa, p = 2.5 GPa, p = 2.75 GPa Hertz peak stress, respectively.

MULTIAXIAL HIGH CYCLE FATIGUE LIMIT CRITERIA
Numerical analysis of the high-cycle RCF process requires knowledge of multiaxial stress-based failure criteria.
In this paper we mainly focus on the multiaxial criteria which use the critical plane concept. In the critical plane concept the orientation of the plane also determines the orientation of the crack. Article [9] defines the critical amounts for pitting, spalling and case crushing. The critical plane is always determined by the maximum shear stress amplitude that occurs in the repetitive cycles, i.e., that plane is the critical plane in which the maximum shear stress amplitude arises during a cycle.
Based on the articles [9]- [11] the Dang Van, Matake and Findley criteria seem to be very promising.
Equivalent stress quantity for the Dang Van criterion [MPa]: where τ(t) is maximal shear stress, MPa; σ h (t) hydrostatic stress, MPa; aDV material parameter; t time. Equivalent stress quantity for the Matake criterion, MPa: where b is material parameter; τ a twice of shear stress amplitude in the critical plane, MPa; σ n, max maximal normal stress in the critical plane, MPa. Equivalent stress quantity for the Findley criterion, MPa: where k is material parameter. The k parameter can be defined by fatigue limits as follows: where σ F , MPa is axial fatigue limit in case of unidirectional (tension-compression) loads with stress ratio R = −1, τ F , MPa tension fatigue limit with stress ratio R = −1.
The above equivalent stress quantities typically all contain material parameters that can only be determined by inverse method, with a sufficiently large number of physical examinations and finite element structural analysis.
However, the results of articles [12] and [13] provide an opportunity to redefine the above equivalent quantities which contain only the fatigue limits obtained as a result of fatigue tests with stress ratio R = −1, and thus their material parameters can be determined without performing an inverse engineering method.
Based on the results of the above-mentioned articles [12] and [13], the criteria for a practically infinite lifetime of a given contact body can be formulated in such a way that the equivalent stresses are comparable to the measurable fatigue limits: Based on the above, the Dang Van criterion can be written as follows: where σ h, max is maximal hydrostatic stress in the critical plane, MPa. The Matake criterion can be written as follows: Technical Gazette 28, 5(2021), 1560-1567 The Findley criterion can be written as follows: The criteria listed can also be rewritten to determine the failure at a given finite life span. As stated above, in the case of the machine elements examined here, a long or ultra-high cycle life must be expected, during which one cannot expect a significant (plastic) deformation from cycle to cycle. Thus, for this typically stress-controlled lifetime estimate, one can use the S-N curve, which is typically given by the Basquin equation [14].
where A, B material parameters, N cycle number of lifespan. Thus, the above fatigue criteria can be defined as follows for fixed lifetimes, based on the results of [12] and [13].
Dang Van criterion for N cycle life: Matake criterion for N cycle life: Findley criterion for N cycle: Since we did not have own measurement results, we tried to decipher the values of σ F , τ F , A and B based on literature data.
The results of the widely accepted Murakami [15] article on σ F , establishes a relationship between the axial (tensile-compressive) fatigue limit and the hardness of the material. This relation for less than 600 HV hardness is as follows: where HV is the Vickers hardness. However, at values higher than 600 HV, saturation is observed; the fatigue limit does not increase further [15].
A similar, well-handled equation for the τ F torsional fatigue limit is not known. Typically, its relationship to the axial fatigue limit is investigated in the literature. According to the evidence of literature [16]- [18], it can be said that in the case of ductile materials there is a significant difference between them. However, the more brittle the material, the closer the torsional and tensioncompression fatigue limits get to each other. According to the literature [14] and [17], in the hardest layers in our case, this ratio is 1.25 and 1 + ν, respectively, where ν is the Poisson's factor, which can vary between 0.25 and 0.33.
Based on the literature [17] a value of −0.09 gives a generally good approximation of B. Thus, B = − 0.09 in our work.
The estimation of the A value is well supported by the article by Roessle and Fatemi [19] which establishes the relationship between Brinell hardness and the coefficient of fatigue life A in the hardness range of our interest. This relationship is shown below:

HB 225
where HB is the Brinell hardness.

COUPLED NUMERICAL MODEL SET UP
In the first step the FEM model of the carburization of cylindrical specimen was set up. The specimen is axisymmetric, so we have simplified the task to a twodimensional model to reduce computational time. The geometry and finite element representation used in the modelling is shown in Fig. 4.

Figure 4 FEM model of carburization
The effect of carbon content and incorporating multiphase transformation models were taken into consideration during thermo-metallurgical and mechanical calculations of the heat treatment cycle.
MSC Marc finite element code was used for numerical calculations. Due to axi-symmetry of the geometry the analysis of the process cylindrical coordinate system was made. Four-Gauss-point axisymmetric rectangular elements were used to model the specimens. Based on expected reduction in carbon content fine to coarse mesh was developed.
For the calculation temperature-dependent, elasticplastic isotropic hardening phase-dependent material model was applied. The thermo-mechanical and metallurgical properties were generated by JMatPro software. The average value of the chemical composition detailed in Tab. 1. was set for the calculation, the austenitization temperature was 920 °C and the initial particle size was 10 microns.
In the simulation, the diffusion was determined according to the Tibbetts proposed relationship [20], while the process-specific weight transfer factor was based on the recommendation of the relevant literature.
The results obtained are illustrated in Figs. 5 to 8. below: Figure 5 Total equivalent plastic stain [-] in the end of the heat treatment process  In the second step the contact model of the roller specimens was set also using the MSC Marc software simulating the state of the rolling contact fatigue. This model was 3-dimensional using periodic boundary condition, symmetry plane and rigid counter-roller as it can be seen in Fig. 9 below. The contact model used the results of the carburization model as initial conditions and material properties. Eight-Gauss-point hexagonal elements were used to model the specimens. The mesh is graded from fine to coarse according to the distance from the chamfer at the edge of the rolling contact surface because the highest mechanical effect is expected there. As it can be seen in Figs. 10 to 12 in the first few cycles there is an elastic shakedown and the extent of the maximum plastic strain is rapidly decreased until the steady state elastic response is reached in the whole body of the specimen.

VALIDATION OF THE COUPLED NUMERICAL MODEL
In order to validate our coupled RCF simulation method and choose the most appropriate failure criterion, the measured results were compared to the calculated outputs. In case of Vickers hardness, the comparison in the symmetry plane of the roller specimen perpendicular to rotating axis is illustrated in Fig. 13.

Figure 13
Comparison of calculated and measured Vickers hardness in the middle plane of the roller specimen the end of the heat treatment process The match is convincing, and one can state that the heat treatment part of the coupled numerical simulation is valid.
The comparison of the contact and failure calculating part of the coupled model with experimental results in case of 10000 N load is illustrated in Figs. 14 to 23 below. In case of the calculated failure criteria, the plotted extent is the left-hand side of the (6)-(8) and (10)- (12) divided by its right-hand side value. Thus, values higher than 1 indicate the failure in the given positions in Fig. 17-23 (and also in Figs. 26 to 28 and Figs. 31 to 33).
In Figs. 18 to 23, there is a comparison between the damage values obtained by endurance limit containing criteria and damage values obtained by lifespan dependent Basquin equation containing criteria. In the latter case the N lifespan was 4519100 according to the performed experiments. As it was expected, the criteria that include the endurance limit always provide higher value with almost the same distribution, because in this case, there is practically infinite time to wait for the failure.
The specimens of experiments loaded by 10000 N show that the spalling and crack propagation occurred at the chamfer in the edge of the rolling contact area (Figs. 14 to 15). This failure was indicated by all the criteria. In Figs. 16 to 17 can be seen that pitting lower than 10 μm appeared on the contact surface. Only the Matake criteria indicate that there is a perceptible increase in the damage value towards the surface, but their values are also much less than 1 at the surface. The difference may be due to the fact that our FEM model cannot take into account the effect of surface roughness and inhomogeneity of surface microstructure. In case of load of 12500 N, we also obtained spalling on chamfer and pitting on contact area (Fig. 24). Similarly to the 10000 N case, the depth of the pitting is less than 10 µm (Fig. 25).
The calculated damage values for only the lifespan driven criteria can be seen in Figs. 26 to 28. Each of them predicted spalling on the edge of the chamfer because their values were higher than 1. Matake criterion showed the higher damage value on the contact surface but values of Findley criterion also have got some increase there. But these values still cannot bridge the gap resulting from surface roughness and inhomogeneity of microstructure. In case of load of 15000 N, we also obtained spalling on chamfer and pitting on contact area (Fig. 29). We also experienced a significant difference when crushing appeared and because of it, the roller specimen was broken as it can be seen in Fig. 30. According to the measurements, our models indicate failure below the carburized layer, i.e. 1.5-2 mm below the surface which can be considered as case crushing. And as before, on the edge of the chamfer, the failure is also indicated by them. Overall, we can say that Matake criterion shows more strongly the damage to places where the failure actually appears.
On the other hand, at a load of 15000 N, Matake criterion already shows failure where it does not exist (300 μm under the surface), and the part of the body which designated for case hardening failure is more extensive in the case of Dang Van (300-400 μm wider). But in the former case, the damage peak coincided with a sudden increase in the calculated hardness gradient. In reality, this place may have shifted towards the surface because of the possible surface grinding and polishing.

CONCLUSION
In this work, we created a coupled numerical simulation method for RCF analysis which takes into account the effects of results of the heat treatment process of the manufacturing of specimen/machine parts, such as hardness, mechanical properties, residual stresses in the carburized layers of the body. Comparing the measured and calculated results, we obtained that the coupled numerical model is valid. The Matake criterion based on the Basquin approach proved to be the best and most conservative predictor of fatigue failure.