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

Abstract

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) K_I and K_II 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 x₀/b = −0.895. Furthermore, the pit shapes and variation of stress intensity factor are determined for various combinations of initial crack length a₀ 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.

Keywords

Lubricated contact pair model micro-pitting contour integral MTS criterion damage characteristics

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³ 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 Glodež’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–9 Ding and Gear¹⁰ 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,¹¹ 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.¹⁴ 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.¹⁵ and Weibring et al.¹⁶ 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–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 β and length a₀, 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:

p (x) = \frac{2 F_{N}}{π b^{2}} \sqrt{b^{2} - x^{2}},

(1)

p_{0} = \sqrt{\frac{F_{N} E^{*}}{2 π R^{*}}},

(2)

q (x) = fp (x),

(3)

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₀ 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.

Figure 1.

(a) Load distribution in the contact area and (b) a crack passing through the oil film and contacting the gear.

The MTS criterion and contour integral method for crack propagation

The MTS criterion is typically used for pitting crack propagation.²⁰ The criteria are as follows: The crack propagation angle θ_C should meet the requirements:

{[K_{I} \sin θ + K_{II} (3 \cos θ - 1)]}_{θ = θ_{C}} = 0,

(4)

where K_I and K_II are mode I and mode II SIFs, respectively.

And the equivalent SIF K_σ in this direction can be determined as:

K_{σ} = \cos \frac{θ}{2} (K_{I} \cos^{2} \frac{θ}{2} - \frac{3}{2} K_{II} \sin θ) .

(5)

If K_σ ≥ 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 Γ 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:

J = \frac{1}{E'} ({K_{I}}^{2} + {K_{II}}^{2}),

(6)

where E′ is defined in terms of E (Young’s modulus) and μ (Poisson’s ratio) of plane strain as:

E' = \frac{E}{1 - μ^{2}} .

(7)

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:

I = \frac{2}{E'} (K_{I} {K_{I}}^{aux} + K_{II} {K_{II}}^{aux}),

(8)

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:

K_{I} = \frac{E'}{2} I_{1}, K_{II} = \frac{E'}{2} I_{2} .

(9)

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 Figure 2 can be established, which can greatly improve the calculation efficiency. The initial crack length and germination angle are a₀ = 20 μm and β = 20°, respectively. The gears are made of flame-hardened AISI 4130 steel with Young’s modulus E = 206 GPa, Poisson’s ratio μ = 0.3, plane strain fracture toughness K_IC = 83.8 MPa m^1/2, and threshold stress intensity factor K_th = 8.5 MPa m^1/2.

Figure 2.

Equivalent contact finite element model that includes oil film thickness and the initial crack.

The equivalent contact model is discretized with quadrilateral-dominated elements. There are singularities in the crack tip region, and other regions are not 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^−1/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.²¹ First, the equivalent radius R^* is equal to 7.1 mm. Then, the Hertzian contact pressure distribution p(x) with a maximum value p₀ = 1550 MPa, and the half-length of the contact area with 0.2 mm, are estimated by using Hertzian contact theory.⁸ ISO-VG-220 lubricating oil, with kinematic viscosity v = 220 mm²/s, density ρ = 900 kg/m³, pressure viscosity coefficient α = 0.018 mm²/N, and entrance velocity u = 5 m/s, is used to calculate the EHL pressure.⁴

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).

Figure 3.

(a) Pressure distribution and (b) different load conditions.

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₀/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 Zafošnik et al.⁸ To verify the correctness of the contact pair model, the Hertzian pressure distribution in Zafošnik et al.⁸ 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 θ_C at the crack tip are calculated by using the contour integral in ABAQUS. Then, the equivalent SIF K_σ is obtained from equation (5). If K_σ ≥ 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 μm (the length of each propagation step is Δa = 10% a₀) every time step in the direction of the propagation angle. When K_σ ≥ 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 μm). Substituting K_I, K_II, and θ_C from the table into equation (4) for the propagation angle θ_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.

Table 1.

Comparison between present results and Zafošnik et al.⁸

	a/μm	K _I/MPa m^1/2	K _II/MPa m^1/2	K _σ/MPa m^1/2	θ_C/°	Equation (4)/10⁻³
Zafošnik et al.⁸	20	8.85	4.18	11.13	−39	4.08
Present results		8.65	4.04	10.84	−38.74	0.46
Error (%)		2.26	3.35	2.61	0.67	/
Zafošnik et al.⁸	22	15.67	3.29	16.64	−22	8.78
Present results		15.50	3.25	16.45	−21.95	0.63
Error (%)		1.08	1.22	1.14	0.23	/
Zafošnik et al.⁸	24	23.19	3.59	23.98	−15	811.01
Present results		26.65	3.64	27.37	−15.01	5.40
Error (%)		14.92	1.39	14.14	0.07	/
Zafošnik et al.⁸	26	53.32	5.66	54.2	−12	136.91
Present results		57.83	6.52	58.91	−12.56	3.91
Error (%)		8.46	15.19	8.69	4.67	/

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.

Table 2.

Numerical results for the moving load considering EHL conditions.

a/μm	Load case	K _I/MPa m^1/2	K _II/MPa m^1/2	Max K_σ/MPa m^1/2	θ_C/°
20	I	0.16	−0.06	10.27 (load case IV)	−33.88 (load case IV)
	II	0.87	0.24
	III	8.31	3.26
	IV	8.74	3.27
	V	8.72	3.05
	VI	8.11	2.55
	VII	7.35	2.05
22	I	0.12	−0.05	14.70 (load case IV)	−19.35 (load case IV)
	II	1.29	0.15
	III	13.46	2.59
	IV	14.05	2.54
	V	13.84	2.28
	VI	12.67	1.80
	VII	11.28	1.32
24	I	0.14	−0.02	22.67 (load case IV)	−14.66 (load case IV)
	II	1.96	0.19
	III	21.23	3.00
	IV	22.10	2.94
	V	21.76	2.66
	VI	19.81	2.10
	VII	17.58	1.56
26	I	0.25	0.01	41.42 (load case IV)	−13.17 (load case IV)
	II	3.57	0.37
	III	38.99	4.84
	IV	40.59	4.81
	V	39.98	4.46
	VI	36.40	3.70
	VII	32.31	2.93
28	I	0.69	0.07	113.73 (load case IV)	−12.40 (load case IV)
	II	9.85	1.05
	III	107.3	12.23
	IV	111.7	12.43
	V	110	11.87
	VI	100.3	10.31
	VII	89.04	8.69

The micro-pitting crack is composed mainly of mode I cracks and some mode II cracks. Simultaneously, the value of K_σ in load case IV is the maximum. When the crack length is up to 28 μm, the crack directly propagates to the tooth surface along the angle of the last propagation step to form micro-pitting because K_σ 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 μm. A contact pair model with h_c = 2 μm 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_σ with crack length a is shown in Figure 4(g). K_σ increases with crack length a in the early stage of crack propagation (a = 20, 22, 24 μm); then, it fluctuates in the middle and late stages (a = 26, 28, 30 μm) due to the obvious increase in pressure inside the crack and the contact effect, thus affecting the change in K_σ.

Figure 4.

Crack propagation process for h_c = 2 μm: (a) a = 20 μm, (b) a = 22 μm, (c) a = 24 μm, (d) a = 26 μm, (e) a = 28 μm, (f) a = 30 μm, and (g) relationship between K_σ and a.

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 μm 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_σ values at the initial crack tip under 12 initial conditions are shown in Table 3.

Figure 5.

Crack propagation in a cross-section of micro pitted gear.

Table 3.

K _σ/MPa m^1/2 values at the initial crack tips for different conditions.

β	a ₀
	10 μm	15 μm	20 μm	25 μm
10°	8.02	11.40	13.57	14.30
20°	3.85	7.01	10.19	12.06
30°	2.80	4.83	6.82	8.62

K _σ increases with initial length a₀ and decreases with a larger germination angle β. 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_σ is larger than 8.5 MPa m^1/2. The results for the six crack propagation cases are further shown in Figure 6.

Figure 6.

(a) Relationship between K_σ and a, (b) relationship between θ_C and a, and (c) comparison of crack propagation paths for different parameters.

According to the above graphs, the following crack propagation laws can be obtained. When β 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₀. In addition, their paths of crack propagations to the edge of the tooth surface are approximately the same. Because a₀ 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 β, 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 micro-pitting 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_σ increases first, then decreases and then increases again. Meanwhile, θ_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_σ and θ_C. In the early stage of crack propagation, the long initial crack withstands much pressure, and the contact area is small. Therefore, K_σ 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_σ 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_σ increases again. Meanwhile, θ_C decreases because the internal pressure of the crack is high enough to reverse the crack propagation direction.

The relationship between K_σ, θ_C, and a can provide theoretical support for contact fatigue life analysis.

Some factors that affect the formation of micro-pitting 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 micro-pitting 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:

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₀/b = −0.895 used in the range of calculations in this paper.

For situations where the initial crack tip K_σ 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.

When the germination angle β is constant, the total length of final crack propagation and micro-pitting area increase with increasing a₀, and the paths of crack propagation are approximately the same.

When the initial length a₀ is constant, the total length of the final crack propagation and micro-pitting area increase with the germination angle β.

If the crack propagates in the deeper area below the gear tooth surface, K_σ may fluctuate significantly, and θ_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.

Footnotes

Handling Editor: James Baldwin

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.

Funding

The author(s) disclosed receipt of the following financial support for the research, authorship, and/or publication of this article: This project is supported by the National Natural Science Foundation of China [Grant No. 51807019, No. 51705057], the China Postdoctoral Science Foundation [Grant No. 2020 M672238], the Chongqing Special Postdoctoral Science Foundation [Grant No. XmT2018030] and the Open Projects of State Key Laboratory for Strength and Vibration of Mechanical Structures [Grant No. SV2020-KF-15].

ORCID iD

Yang Xiao

References

Höhn

Michaelis

. Influence of oil temperature on gear failures. Tribol Int 2004; 37: 103–109.

Chen

Liu

. Study of the influence of lubrication parameters on gear lubrication properties and efficiency. Ind Lubr Tribol 2016; 68: 647–657.

Aslantas

Tasgetiren

. A study of spur gear pitting formation and life prediction. Wear 2004; 257: 1167–1175.

Fajdiga

Flašker

Glodež

, et al. Numerical modelling of micro-pitting of gear teeth flanks. Fatigue Fract Eng Mater Struct 2003; 26: 1135–1143.

Fajdiga

Flašker

Glodež

. The influence of different parameters on surface pitting of contacting mechanical elements. Eng Fract Mech 2004; 71: 747–758.

Glodež

Aberšek

Flašker

, et al. Evaluation of the service life of gears in regard to surface pitting. Eng Fract Mech 2004; 71: 429–438.

Fajdiga

Glodež

Kramar

. Pitting formation due to surface and subsurface initiated fatigue crack growth in contacting mechanical elements. Wear 2007; 262: 1217–1224.

Zafošnik

Glodež

Ulbin

, et al. A fracture mechanics model for the analysis of micro-pitting in regard to lubricated rolling–sliding contact problems. Int J Fatigue 2007; 29: 1950–1958.

Glodež

Potočnik

Flašker

, et al. Numerical modelling of crack path in the lubricated rolling–sliding contact problems. Eng Fract Mech 2008; 75: 880–891.

10.

Ding

Gear

. Spalling depth prediction model. Wear 2009; 267: 1181–1190.

11.

Anisetti

. A tribo-dynamic contact fatigue model for spur gear pairs. Int J Fatigue 2017; 98: 81–91.

12.

Cao

Wang

, et al. Effect of contact path on the mixed lubrication performance, friction and contact fatigue in spiral bevel gears. Tribol Int 2018; 123: 359–371.

13.

Wang

Zhang

, et al. Transient behaviors of friction, temperature and fatigue in different contact trajectories for spiral bevel gears. Tribol Int 2019; 141: 105965.

14.

Liu

Zhu

, et al. Study on the gear fatigue behavior considering the effect of residual stress based on the continuum damage approach. Eng Fail Anal 2019; 104: 531–544.

15.

Lai

Lohmann

, et al. A model to predict initiation and propagation of micro-pitting on tooth flanks of spur gears. Int J Fatigue 2019; 122: 106–115.

16.

Weibring

Gondecki

Tenberge

. Simulation of fatigue failure on tooth flanks in consideration of pitting initiation and growth. Tribol Int 2019; 131: 299–307.

17.

Chen

Cao

. Research on dynamics and fault mechanism of spur gear pair with spalling defect. J Sound Vib 2012; 331: 2097–2109.

18.

Morales-Espejel

Gabelli

. A model for gear life with surface and subsurface survival: tribological effects. Wear 2018; 404–405:133–142.

19.

Meng

Wang

Shi

. A novel evolution model of pitting failure and effect on time-varying meshing stiffness of spur gears. Eng Fail Anal 2020; 120: 105068.

20.

Erdogan

Sih

. On the crack extension in plates under plane loading and transverse shear. J Basic Eng 1963; 85: 519–525.

21.

Huang

. Numerical calculation methods of elastohydrodynamic lubrication. Beijing: Tsinghua University Press, 2013, pp. 64–83.

22.

Evans

Snidle

Sharif

, et al. Analysis of micro-elastohydrodynamic lubrication and prediction of surface fatigue damage in micropitting tests on helical gears. J Tribol 2013; 135: 011501.