A lubrication contact pair model for simulating gear micro-pitting damage characteristics based on contour integral

A new two-dimensional finite element model of a lubricated contact pair, based on a contour integral, is proposed to investigate the formation of micro-pitting on gear tooth surfaces. Meanwhile, the contact properties and elasto-hydrodynamic lubrication (EHL) conditions of the gears are considered in the lubricated contact pair model. Then, the stress intensity factors (SIFs) KI and KII and the propagation angle θ C at the crack tip are analyzed by ABAQUS software. Next, the equivalent SIF Kσ can be calculated according to the maximum tangential stress (MTS) criterion, which is often used as the criterion for crack propagation. Considering the effect of a moving contact, the crack more easily propagates under the load x0/b = −0.895. Furthermore, the pit shapes and variation of stress intensity factor are determined for various combinations of initial crack length a0 and angle β. The results show that longer germinated cracks propagate in areas that are deeper below the tooth surface. And the total length of final crack increases with the initial length and germination angle. These research results provide theoretical support for contact fatigue life analysis and meshing stiffness calculations of micro-pitting gears.


Introduction
Shortly after a gear is put into service, many micropits may form on the surfaces of gear teeth and cause gear wear. As a result, the parameters of gear tooth surfaces and dynamic characteristics between tooth surfaces change, which leads to failure of the whole mechanism. Therefore, it is necessary to study the formation of the micro-pitting starting from crack initiation and propagation to the surface and reveal the influence of EHL conditions on the damage characteristics of micro-pitting, which provides a theoretical reference for slowing gear wear and prolonging gear lifetime.
Experimental studies on gear micro-pitting indicate that the formation of micro-pitting is affected by many factors, such as gear materials, surface impurities, roughness, loads, and lubricating oil temperatures, 1,2 which makes it difficult to build a numerical model by taking these factors into account. Current research mainly focuses on the equivalent contact fatigue model. Aslantas and Tasgetiren 3 proposed a pitting failure model based on finite element analysis and linear elastic fracture mechanics, and they used the MTS criterion to determine the crack propagation direction. An equivalent computational model of gear tooth surface contact fatigue was established by Glodezˇ's team. They predicted crack propagation paths by virtual crack growth and contour integral methods, and the crack propagation direction was determined by the MTS criterion and modified MTS criterion. Finally, the relationship between the SIFs and crack length a in the formation of micro-pitting was obtained; furthermore, the number of stress cycles required for a crack to propagate from initiation to critical crack length was calculated. In addition, the formation of micro-pitting for different loading conditions, initial cracks, and mechanical properties was studied by an equivalent computational model. [4][5][6][7][8][9] Ding and Gear 10 proposed a predictive model of gear spalling depth and used the model to predict the spalling depth of spur gears with different geometries and compositions.
In recent years, scholars have proposed new models to study pitting of gears. A dynamic contact fatigue model under high speed conditions was established by Li and Anisetti, 11 and the effects of torque, surface roughness, and lubrication temperature on contact fatigue of gears were studied by using the model. Researchers established a hybrid EHL model and used it to study the effects of lubrication characteristics, friction coefficient, temperature, and contact fatigue on the lifetime of spiral bevel gears. 12,13 He et al. 14 used a numerical model of gear fatigue damage that considered the initial residual stress value of 380 MPa to simulate the contact fatigue of a gear tooth surface. Xu et al. 15 and Weibring et al. 16 proposed a numerical model for predicting the initiation and propagation of micro-pitting on the tooth surface of spur gears.
Previous studies on crack propagation and pitting of gears mainly based on Hertz theory and EHL theory to establish a numerical equivalent contact model of a gear tooth surface with micro-cracks. However, the influence of the contact between a crack and the opposite tooth surface was not considered during crack propagation. In addition, it is necessary to clarify the change in tooth surface morphology on the scale characteristics of micro-pitting to establish an analytical model of the influence of micro pitting, whether to study of normal contact stiffness or predict the contact fatigue life of tooth surfaces. [17][18][19] Unlike other studies, this paper considers the influence of the EHL condition of the gear as well as the contact between the crack and the opposite tooth surface on the formation of micro-pitting. Based on a contour integral, a numerical finite element model of lubricated contact pairs is established to study the formation of micro-pitting on gear tooth surfaces. Then, the damage characteristics of micro-pitting, which are influenced by initial cracks with different values of the germination angle b and length a 0 , are obtained. The research results can provide theoretical support for contact fatigue life analysis and meshing stiffness calculation of micro-pitting on gears.

Factors influencing crack propagation in lubricated contacts
According to Hertzian theory, which is widely used to model unlubricated surfaces, the contact geometry of gear tooth surfaces can be transformed into a pair of equivalent contact cylinders. The distribution of the normal contact pressure p(x) and tangential contact load q(x) can be analytically determined by: where F N is the normal force per unit width of equivalent cylinders, b is the half-length of the contact area, f is the friction coefficient, E * is the equivalent Young's modulus, R * is the equivalent radius, and p 0 is the maximum contact pressure. However, in a lubricated contact, the contact pressure between gear pairs differs, as shown in Figure 1(a), which is most apparent in the inlet area, outlet area, and pressure spike. At the same time, the gear pair is separated by an oil film under EHL conditions, which can transfer the contact pressure. Therefore, the normal load between gear pairs is equal to the value of the EHL pressure. In addition, a certain pressure at the crack mouth forms inside the crack and forces the crack to propagate in the forward direction. The crack mouth area warps and deforms during crack propagation. When the deformation is large, the crack mouth passes through the oil film and contacts the opposite gear, as shown in Figure 1(b). The contact may affect the characteristics of crack propagation.

The MTS criterion and contour integral method for crack propagation
The MTS criterion is typically used for pitting crack propagation. 20 The criteria are as follows: The crack propagation angle u C should meet the requirements: where K I and K II are mode I and mode II SIFs, respectively. And the equivalent SIF K s in this direction can be determined as: If K s ø K IC , then the brittle fracture of the part occurs due to one-time loading. However, in the process of fatigue crack propagation, the critical SIF K th replaces K IC as the criterion of fatigue crack propagation. K IC and K th are the plane strain fracture toughness and threshold stress intensity factor, respectively. Similar to Young's modulus E, they are specific performance parameters of materials.
In addition, the stress, strain, and displacement in a crack tip region are integrated along a closed-loop G to obtain a constant J-integral independent of the path, which is called the contour integral method. The relationship between the J-integral and stress intensity factor is as follows: where E# is defined in terms of E (Young's modulus) and m (Poisson's ratio) of plane strain as: For mixed-mode cracks, a hypothetical auxiliary field at the crack tip, combined with the interaction integral method, is used to calculate the SIFs. The relationship between I and the SIFs of various fields can be obtained by: where K I aux and K II aux are mode I and II SIFs of the auxiliary field. K I aux = 1, K II aux = 0 and K I aux = 0, K II aux = 1 are introduced into the above formula in turn. The mode I and II SIFs can be obtained through the calculation of two interaction integrals: Simulation of the micro-pitting of gears caused by a lubricated contact

Numerical model
Based on the contour integral and MTS criterion, the propagation of gear pitting cracks is numerically simulated by means of ABAQUS software of French Dassault Company. The radius on the meshing line changes with the meshing position. In addition, the contact fatigue strength of the tooth surface near the pitch line is usually the weakest, so the area near the pitch line was selected to study micro-pitting. According to Hertzian contact theory, the meshing of gears can be transformed into two cylindrical contacts, and then it can be further simplified as a plane strain problem according to elastic mechanics. Due to the symmetry principle, the circular section can be simplified as a semicircle, and the contact gap between the two semicircles corresponds to the oil film thickness. Therefore, the two-dimensional finite element model in  affected, so 36 collapsed quadrilateral quarter-point finite elements are used around the crack tip, which can not only meet the accuracy requirements but also simulate the r 21/2 stress singularity at the crack tip. The other regions of the semicircle are discretized with sparse elements. The constraint in the y-direction is applied to the whole bottom of the half-cylinder, and the constraint in the x-direction is applied at the center of the half-cylinder at the same time. The boundary conditions and the initial crack are shown in Figure 2. The corresponding normal and tangential loads are applied to the contact area of the gear pair, and the friction factor affecting the tangential load is f = 0.04. The friction coefficient of the part where the pitting crack passes through the oil film and contacts the opposite tooth surface is taken as 0.1.

Contact pressure distributions
The EHL pressure of the gear pair is calculated by the numerical method of isothermal line contact elasto-hydrodynamic lubrication. 21 First, the equivalent radius R * is equal to 7.1 mm. Then, the Hertzian contact pressure distribution p(x) with a maximum value p 0 = 1550 MPa, and the half-length of the contact area with 0.2 mm, are estimated by using Hertzian contact theory. 8 ISO-VG-220 lubricating oil, with kinematic viscosity v = 220 mm 2 /s, density r = 900 kg/m 3 , pressure viscosity coefficient a = 0.018 mm 2 /N, and entrance velocity u = 5 m/s, is used to calculate the EHL pressure. 4 When solving the EHL pressure distribution, the Hertzian contact pressure distribution should be taken as the initial pressure distribution to calculate the film thickness and lubricant viscosity, and the new pressure distribution can be solved by the Reynolds equation. Then, the previous pressure distribution is iteratively corrected, and the elastic deformation and film thickness are recalculated. The iteration is repeated until the pressure difference between the two iterations is very small. Finally, the calculated EHL pressure distribution, film thickness, and corresponding Hertzian pressure distribution are shown in Figure 3(a).
Previous studies have shown that cracks propagate most easily when the moving loaded contact region reaches the crack mouth. 8,9 Therefore, seven loading configurations were considered to investigate the effect of a moving contact on the stress field around the crack tip, as shown in Figure 3(b). Each loading configuration has the same normal p(x) and tangential q(x) contact loading distributions but acts at different positions with respect to the initial crack. x 0 /b shows the distance from the crack mouth to the position of the maximum contact pressure.

Validation of the model and determination of the load conditions
Validation of the model. The elasto-hydrodynamic lubrication pressure and the influence of the contact  between the crack and the opposite tooth surface were not considered in Zafosˇnik et al. 8 To verify the correctness of the contact pair model, the Hertzian pressure distribution in Zafosˇnik et al. 8 is applied to the model. Furthermore, it does not contact the opposite gear tooth surface after crack deformation; that is, the oil film thickness between the contact pair models in this calculation is large enough. The K I , K II , and propagation angle u C at the crack tip are calculated by using the contour integral in ABAQUS. Then, the equivalent SIF K s is obtained from equation (5). If K s ø K th , the crack extends with the propagation angle. To obtain the crack propagation path accurately, the length of each propagation should not be too long. Therefore, the length of the crack increases by 2 mm (the length of each propagation step is Da = 10% a 0 ) every time step in the direction of the propagation angle. When K s ø K IC , the crack propagates directly to the tooth surface along the propagation angle, and spalling occurs.
In the reported references, the effect of crack propagation (not yet extending to the tooth surface to form pits) on the EHL pressure is not quantitatively analyzed. In addition, because the crack deformation area is very small relative to the whole contact area, it is assumed that the crack deformation has little effect on the contact pressure distribution of the whole contact surface.
The comparison between the calculated results and the literature is shown in Table 1. The two results are in good agreement, and the error is within 4% at the early stage of crack propagation, while the error increases when the crack is about to propagate to the tooth surface (a = 24, 26 mm). Substituting K I , K II , and u C from the table into equation (4) for the propagation angle u C , the results obtained in this paper are more in line with the requirements of the crack propagation angle. Thus, the contact pair model and method in this paper are correct and can be further applied to crack propagation calculations under EHL.
Determination of EHL load conditions. The influence of the EHL pressure distribution is investigated with the above numerical model. After the computation with seven conditions in Figure 3(b), the results are shown in Table 2.
The micro-pitting crack is composed mainly of mode I cracks and some mode II cracks. Simultaneously, the value of K s in load case IV is the maximum. When the crack length is up to 28 mm, the crack directly propagates to the tooth surface along the angle of the last propagation step to form micro-pitting because K s is larger than the plane strain fracture toughness K IC . Therefore, subsequent crack propagation parametric studies are carried out under the load distribution of load case IV.

Parametric study of crack propagation
The crack may contact the opposite gear tooth surface during propagation, so it is necessary to consider the situation in which the crack passes through the oil film and contacts the opposite tooth surface. As shown in Table 2, the crack propagates more easily under load case IV. According to Figure 3(a), the oil film thickness in the crack propagation area is calculated as h c ' 2 mm. A contact pair model with h c = 2 mm is established, and crack propagation is shown in Figure 4(a) to (f). The results show that the micro-pitting crack contacts the opposite tooth surface in the early stage of propagation, and the lubricated contact model is closer to the real situation. The variation in K s with crack length a is shown in Figure 4(g). K s increases with crack length a in the early stage of crack propagation (a = 20, 22, 24 mm); then, it fluctuates in the middle and late stages (a = 26, 28, 30 mm) due to the obvious increase in pressure inside the crack and the contact effect, thus affecting the change in K s .
In this scanning electron microscopic image (SEM), as shown in Figure 5, micro-pitting is characterized by shallow pits with depths that usually do not exceed 10 mm and propagation angles that range from 10°to 30°from the tooth surface. 16,22 Parametric studies are performed by means of numerical simulations to determine the initial conditions that influence crack propagation. The K s values at the initial crack tip under 12 initial conditions are shown in Table 3. K s increases with initial length a 0 and decreases with a larger germination angle b. That is, if the length of the germinated crack on the tooth surface is short, it propagates near the edge of the tooth surface, and a germinated crack with a longer length propagates in the deeper area below the tooth surface, which is consistent with crack propagation shown in Figure 5.
Meanwhile, six initial cracks propagate because the value of K s is larger than 8.5 MPa m 1/2 . The results for the six crack propagation cases are further shown in Figure 6.
According to the above graphs, the following crack propagation laws can be obtained. When b is constant (cases 1-3 or cases 4 and 5), the total length of final crack propagation and the micro-pitting area increase with increasing a 0 . In addition, their paths of crack propagations to the edge of the tooth surface are approximately the same. Because a 0 is a constant (cases 2 and 4 or cases 3, 5, and 6), the total length of the final crack propagation increases with the germination angle b, and the micro-pitting area is also larger. Meanwhile, the crack closer to the edge of the tooth surface can propagate forward more easily and may form micropitting pits under fewer stress cycles. If the crack propagates in the deeper area of the gear tooth surface, such as in cases 5 and 6, K s increases first, then decreases and then increases again. Meanwhile, u C increases first and then decreases. This is because the initial length crack induces the internal pressure and the contact effect to affect the values of K s and u C . In the early stage of crack propagation, the long initial crack withstands much pressure, and the contact area is small. Therefore, K s can increase during this period. In the middle stage of crack propagation, the increasing contact area counters the effect of the internal pressure of the crack. As a result, K s tends to decrease. Finally, the crack tip is closer to the edge of the surface with further propagation, which can make the stress singularity at the crack tip more obvious, so K s increases again. Meanwhile, u C decreases because the internal pressure of the crack is high enough to reverse the crack propagation direction.
The relationship between K s , u C , and a can provide theoretical support for contact fatigue life analysis.
Some factors that affect the formation of micropitting are ignored in the contact pair model, such as surface impurities, roughness, and the influence of crack deformation on the contact pressure. By comparing the micro-pitting pits in Figures 5 and 6(c), it can be seen from the results that the shape and size of micropitting pits determined by numerical simulation are close to the actual situation, which not only shows that these neglected factors have little influence on the present results, but also proves the validity of the model and method.

Conclusions
In this paper, a new two-dimensional finite element numerical model of lubricated contact pairs, based on contour integrals, is proposed. Hertzian theory, contact surface friction, and EHL conditions are all considered in the model. The main conclusions are as follows: (1) The contact pair model with an EHL pressure distribution is more consistent with the actual operating conditions of a gear, and the crack propagates more easily for the load x 0 / b = 20.895 used in the range of calculations in this paper. (2) For situations where the initial crack tip K s exceeds K th , if the germinated crack is short, it propagates near the edge of the tooth surface, and longer germinated cracks propagate in areas that are deeper below the tooth surface. (3) When the germination angle b is constant, the total length of final crack propagation and micro-pitting area increase with increasing a 0 , and the paths of crack propagation are approximately the same. (4) When the initial length a 0 is constant, the total length of the final crack propagation and micro-pitting area increase with the germination angle b. (5) If the crack propagates in the deeper area below the gear tooth surface, K s may fluctuate significantly, and u C first increases and then decreases.
The micro-pitting pit morphology and corresponding mechanical properties can be obtained quickly and  accurately by the lubricated contact pair model, which provides a theoretical reference for fatigue lifetime analysis and meshing stiffness calculations of micro-pitting gears. The next step is to focus on crack growth behaviors under various roughnesses and temperatures.

Declaration of conflicting interests
The author(s) declared no potential conflicts of interest with respect to the research, authorship, and/or publication of this article.