A method and system for predicting the contact fatigue life of spiral bevel gears

Abstract

The present invention discloses a method and system for predicting the contact fatigue life of spiral bevel gears, including calculating equivalent meshing characteristic parameters according to blank parameters and machining parameters; calculating first characteristic parameters of the load-bearing contact points on the tooth surface of the bevel gear according to the equivalent meshing characteristic parameters; performing numerical calculation through the non-Newtonian fluid thermo-elastohydrodynamic lubrication equation based on the first characteristic parameters and the true rough surface topography to obtain the normal pressure in the contact area; performing numerical calculation through the non-Newtonian fluid thermo-elastohydrodynamic lubrication equation based on the first characteristic parameters and the true rough surface topography to obtain the tooth surface friction force; calculating the subsurface stress-strain field according to the normal pressure in the contact area, the residual stress distributed along the tooth surface depth direction, the microhardness, and the tooth surface friction force; calculating the contact fatigue life of the spiral bevel gear according to the subsurface stress-strain field, providing a theoretical method and data support for the anti-fatigue design and manufacturing of spiral bevel gears, and improving the accuracy of predicting the contact fatigue life of spiral bevel gears.

Description

Technical Field

[0001] The present invention relates to the technical field related to the prediction of the contact fatigue life of spiral bevel gears, and particularly to a method and system for predicting the contact fatigue life of spiral bevel gears. Background Art

[0002] Due to the advantages of strong load-bearing capacity, high contact ratio and smooth transmission, spiral bevel gears have been widely used in high-speed and heavy-load working conditions in fields such as aerospace, automobiles and ships. Due to the complex contact geometry, heavy-load working conditions and tooth surface roughness of spiral bevel gear pairs, the actual lubrication state of the gears is a harsh mixed lubrication state, and the external load is jointly borne by the fluid lubrication area and the micro-convex body contact area. Too low lubricating oil film will cause the direct contact of the tooth surface micro-convex bodies, increasing the stress concentration on the tooth surface; and it is not easy to form a lubricating oil film due to the increase in tooth surface temperature, exacerbating the degree of gear wear. Therefore, studying the mixed lubrication performance of bevel gears has important engineering significance for improving the service life of gears.

[0003] Contact analysis and lubrication analysis are the basis for predicting the contact fatigue life of gears. Due to the complex transmission form and heavy-load working conditions of spiral bevel gears, the change of the actual contact trajectory during loading has a non-negligible impact on the gear lubrication performance. The fatigue life prediction of spiral bevel gears involves multi-disciplinary coupling such as interface contact mechanics, hydrodynamic lubrication mechanics, and fatigue mechanics, and the calculation amount is huge, bringing great challenges to the establishment of a comprehensive and systematic method for predicting the contact fatigue life of spiral bevel gears. On the other hand, the traditional gear life assessment method is based on meeting strength and life theories, lacking a deep understanding of the gear fatigue failure mechanism and the surface integrity action mechanism. The gear surface integrity parameters are crucial for the contact fatigue performance of gears. Therefore, it is urgent to consider the surface integrity and actual load-bearing contact characteristics of spiral bevel gears. Summary of the Invention

[0004] The present invention aims to at least solve the technical problems existing in the prior art. For this reason, the present invention proposes a method and system for predicting the contact fatigue life of spiral bevel gears, comprehensively considering the influence of tooth surface topography-lubrication coupling, residual stress and microhardness on the gear fatigue performance, providing a theoretical method and data support for the anti-fatigue design and manufacture of spiral bevel gears, and improving the accuracy of predicting the contact fatigue life of spiral bevel gears.

[0005] In the first aspect of the present invention, a method for predicting the contact fatigue life of spiral bevel gears is provided, including the following steps:

[0006] Obtain the blank parameters, machining parameters, true rough surface topography, residual stress and microhardness distributed along the tooth surface depth direction of the spiral bevel gear;

[0007] Calculate the equivalent meshing characteristic parameters according to the blank parameters and the machining parameters;

[0008] Calculate the first characteristic parameter of the load-bearing contact point on the tooth surface of the bevel gear according to the equivalent meshing characteristic parameters;

[0009] Perform numerical calculation according to the first characteristic parameter and the true rough surface topography through the non-Newtonian fluid thermo-elastohydrodynamic lubrication equation to obtain the normal pressure in the contact area;

[0010] Perform numerical calculation according to the first characteristic parameter and the true rough surface topography through the non-Newtonian fluid thermo-elastohydrodynamic lubrication equation to obtain the tooth surface friction force;

[0011] Calculate the subsurface stress-strain field according to the normal pressure in the contact area, the residual stress distributed along the tooth surface depth direction, the microhardness, and the tooth surface friction force;

[0012] Calculate the contact fatigue life of the spiral bevel gear according to the subsurface stress-strain field.

[0013] According to the control method of the embodiment of the present invention, it has at least the following beneficial effects:

[0014] This method obtains the blank parameters, machining parameters, true rough surface topography, residual stress distributed along the tooth surface depth direction, and microhardness of the spiral bevel gear; calculates the equivalent meshing characteristic parameters according to the blank parameters and machining parameters; calculates the first characteristic parameter of the load-bearing contact point on the tooth surface of the bevel gear according to the equivalent meshing characteristic parameters; performs numerical calculation according to the first characteristic parameter and the true rough surface topography through the non-Newtonian fluid thermo-elastohydrodynamic lubrication equation to obtain the normal pressure in the contact area; performs numerical calculation according to the first characteristic parameter and the true rough surface topography through the non-Newtonian fluid thermo-elastohydrodynamic lubrication equation to obtain the tooth surface friction force; calculates the subsurface stress-strain field according to the normal pressure in the contact area, the residual stress distributed along the tooth surface depth direction, the microhardness, and the tooth surface friction force; calculates the contact fatigue life of the spiral bevel gear according to the subsurface stress-strain field, comprehensively considers the influence of tooth surface topography-lubrication coupling, residual stress, and microhardness on the fatigue performance of the gear, provides a theoretical method and data support for the anti-fatigue design and manufacturing of the spiral bevel gear, and improves the prediction accuracy of the contact fatigue life of the spiral bevel gear.

[0015] According to some embodiments of the present invention, the equivalent meshing characteristic parameters include equivalent meshing force, equivalent meshing torque, and equivalent meshing point position, and the calculating the equivalent meshing characteristic parameters according to the blank parameters and the machining parameters includes:

[0016] Establish a finite element model according to the blank parameters and the machining parameters;

[0017] Perform a loading contact analysis based on the finite element model to obtain the equivalent meshing force, the equivalent meshing moment, and the position of the equivalent meshing point. Among them, the calculation formula for obtaining the equivalent meshing force according to the loading contact analysis based on the finite element model is as follows:

[0018]

[0019] Among them, F mf is the equivalent meshing force, N is the number of tooth pairs that are simultaneously meshing at any position of the gear, x is the abscissa of the gear, y is the ordinate of the gear, z is the vertical coordinate of the gear, F mx is the transverse equivalent meshing force, F my is the longitudinal equivalent meshing force, F mz is the vertical equivalent meshing force, F mj，i is the coordinate equivalent meshing force of the i-th tooth pair, F mj is the coordinate equivalent meshing force, and j is the sub-coordinate of the gear;

[0020] The calculation formula for obtaining the equivalent meshing moment according to the loading contact analysis based on the finite element model is as follows:

[0021] n j = F mj / F mf (j = x, y, z)

[0022]

[0023] Among them, n j is the direction vector of the equivalent meshing force, M mj is the equivalent meshing moment, and M mj,i is the equivalent meshing moment of the i-th tooth pair;

[0024] The calculation formula for obtaining the position of the equivalent meshing point according to the loading contact analysis based on the finite element model is as follows:

[0025]

[0026] Among them, x m is the abscissa of the position of the equivalent meshing point, y m is the ordinate of the position of the equivalent meshing point, z m is the vertical coordinate of the position of the equivalent meshing point, and the position of the equivalent meshing point is r m (x m , y m , z m ), M mx is the transverse equivalent meshing moment, M my is the longitudinal equivalent meshing moment, and x m，iis the abscissa of the equivalent meshing point position of the i-th tooth pair.

[0027] According to some embodiments of the present invention, the first characteristic parameters include the radius of curvature of the major axis of the contact ellipse, the radius of curvature in the minor axis direction, the relative sliding speed, the entrainment speed and the slide-roll ratio, and the included angle between the entrainment speed and the minor axis of the contact ellipse. Calculating the first characteristic parameters of the load-bearing contact point on the tooth surface of the bevel gear according to the equivalent meshing characteristic parameters includes:

[0028] Performing curvature calculation according to the equivalent meshing point position, the first fundamental form and the second fundamental form of the tooth surface to obtain the principal curvatures (k g1 , k g2 , k p1 , k p2 ) of the tooth surfaces of the large and small gears at the load-bearing contact point and the corresponding principal directions (t g1 , t g2 , t p1 , t p2 );

[0029] Calculating the normal curvature k g1 , k g2 , k p1 , k p2 ) of the large gear and the normal curvature k g1 , k g2 , k p1 , k p2 ) of the small gear and the induced normal curvature k gx , k px , k rx of the load-bearing contact point on the tooth surface of the bevel gear according to the meshing characteristic parameters, the principal curvatures (k g1 , k g2 , k p1 , k p2 ) and the corresponding principal directions (t g1 , t g2 , t p1 , t p2 ). The calculation formulas for the normal curvature of the large gear, the normal curvature of the small gear and the induced normal curvature in the x direction of the load-bearing contact point on the tooth surface of the bevel gear are:

[0030] δ = arccos(t g1 ·t p1 / (|t g1 ||t p1 |))

[0031]

[0032] Where is the angle between any tangent line x and t on the tangent plane of the load-bearing contact point on the tooth surface of the bevel gear, δ is the angle between the t direction of the large gear and the t direction of the small gear, k is the normal curvature of the large gear in the x direction, k is the normal curvature of the small gear, and k is the induced normal curvature; g1 g1 direction and the t direction of the small gear p1 gx is the normal curvature of the large gear in the x direction, k px is the normal curvature of the small gear, k rx is the induced normal curvature;

[0033] Take the maximum and minimum values of the induced normal curvature to obtain the angle corresponding to the maximum value of the induced normal curvature and the angle corresponding to the minimum value of the induced normal curvature where, is the angle between the short axis x-axis of the contact ellipse and t g1 is the angle between the long axis y-axis of the contact ellipse and t g1

[0034] Calculate the length of the short semi-axis and the length of the long semi-axis of the contact ellipse according to the equivalent meshing characteristic parameters and the principal curvatures. Among them, the calculation formulas for calculating the length of the short semi-axis and the length of the long semi-axis of the contact ellipse according to the equivalent meshing characteristic parameters and the principal curvatures are:

[0035]

[0036] where, a is the length of the short semi-axis of the contact ellipse, b is the length of the short semi-axis of the contact ellipse, k a is the first coefficient of the elliptic integral function, k b is the second coefficient of the elliptic integral function, F n is the normal load, A is the first constant related to the shape of the object, B is the second constant related to the shape of the object, E1 is the elastic modulus of the large gear material, E2 is the elastic modulus of the small gear material, μ1 is the Poisson's ratio of the large gear material, and μ2 is the Poisson's ratio of the small gear material;

[0037] Obtain the moving speeds at the load-bearing contact points of the large and small gears according to the basic kinematics of gear transmission. Among them, the calculation formulas for obtaining the moving speeds at the load-bearing contact points of the large and small gears according to the basic kinematics of gear transmission are:

[0038] v p = ω p × r p , v g = ω g × r g

[0039] where, ω g is the angular velocity vector of the large gear, ω p is the angular velocity vector of the small gear, r​​​​g is the radius vector of the large gear, r p is the radius vector of the small gear, v g is the velocity vector of the large gear, v p is the velocity vector of the small gear;

[0040] Calculate the entrainment velocity along the minor axis and major axis of the contact ellipse according to the velocity vector of the large gear and the velocity vector of the small gear, wherein, the calculation formula for calculating the entrainment velocity along the minor axis and major axis of the contact ellipse according to the velocity vector of the large gear and the velocity vector of the small gear is:

[0041]

[0042]

[0043] Wherein, is the tangential plane component of the large gear velocity vector along the minor axis direction of the contact ellipse, is the tangential plane component of the small gear velocity vector along the minor axis direction of the contact ellipse, is the component of the large gear velocity vector along the tooth surface tangential plane, is the component of the small gear velocity vector along the tooth surface tangential plane, is the component of the entrainment velocity along the minor axis direction of the contact ellipse, is the tangential plane component of the entrainment velocity along the major axis direction of the contact ellipse, U e is the entrainment velocity vector, θ e is the angle between the entrainment velocity vector and the minor axis x-axis of the contact ellipse;

[0044] Calculate the entrainment velocity, the angle between the entrainment velocity vector and the minor axis x-axis of the contact ellipse, the sliding velocity, and the angle between the sliding velocity vector and the minor axis x-axis of the contact ellipse according to the entrainment velocity along the minor axis and major axis of the contact ellipse, wherein, the calculation formula for calculating the entrainment velocity, the angle between the entrainment velocity vector and the minor axis x-axis of the contact ellipse, the sliding velocity, and the angle between the sliding velocity vector and the minor axis x-axis of the contact ellipse according to the entrainment velocity along the minor axis and major axis of the contact ellipse is:

[0045]

[0046]

[0047] Wherein, U s is the sliding velocity vector, θ s is the angle between the sliding velocity vector and the minor axis x-axis of the contact ellipse;

[0048] Calculate the slip-roll ratio based on the entrainment speed and the sliding speed, where the calculation formula for calculating the slip-roll ratio based on the entrainment speed and the sliding speed is:

[0049] AKC = U s / U e

[0050] where AKC is the slip-roll ratio of the contact point.

[0051] According to some embodiments of the present invention, the numerical calculation of the normal pressure in the contact area through the non-Newtonian fluid thermo-elastohydrodynamic lubrication equation based on the first characteristic parameter and the true rough surface topography includes:

[0052] Calculate the Reynolds equation for the non-Newtonian fluid thermo-elastohydrodynamic lubrication equation by the separated flow method based on the first characteristic parameter, where the calculation formula for calculating the Reynolds equation for the non-Newtonian fluid thermo-elastohydrodynamic lubrication equation by the separated flow method based on the first characteristic parameter is:

[0053]

[0054] where p is the pressure, h is the oil film thickness, ρ is the lubricating oil density, η is the lubricating oil viscosity, x is the coordinate in the movement direction, y is the coordinate perpendicular to the movement direction, φ x is the dimensionless flow factor one, φ y is the dimensionless flow factor two;

[0055] Calculate the film thickness equation for the non-Newtonian fluid thermo-elastohydrodynamic lubrication equation by the separated flow method based on the first characteristic parameter, where the calculation formula for calculating the film thickness equation for the non-Newtonian fluid thermo-elastohydrodynamic lubrication equation by the separated flow method based on the first characteristic parameter is:

[0056]

[0057] where h0 is the central film thickness when not deformed, R x is the radius of curvature in the x direction, R y is the radius of curvature in the y direction, E′ is the comprehensive elastic modulus, x′ is the additional coordinate corresponding to x, y′ is the additional coordinate corresponding to y, Ω is the solution domain, and S(x, y) is the roughness height value of the tooth surface;

[0058] Calculate the viscosity equation for the non-Newtonian fluid thermo-elastohydrodynamic lubrication equation by the separated flow method based on the first characteristic parameter, where the calculation formula for calculating the viscosity equation for the non-Newtonian fluid thermo-elastohydrodynamic lubrication equation by the separated flow method based on the first characteristic parameter is:

[0059] s0 = β T(T0 - 138) / (lnη0 + 9.67)

[0060]

[0061] Among them, η0 is the ambient viscosity, z0 is the dimensionless viscosity - pressure index, T is the lubricating oil temperature, T0 is the initial temperature of the lubricating oil, s0 is the dimensionless viscosity - temperature index, β T is the viscosity - temperature coefficient, and η is the viscosity;

[0062] Calculate the density equation for the non - Newtonian fluid thermo - elastohydrodynamic lubrication equation by the separated flow method according to the first characteristic parameter. Among them, the calculation formula for calculating the density equation for the non - Newtonian fluid thermo - elastohydrodynamic lubrication equation by the separated flow method according to the first characteristic parameter is:

[0063]

[0064] Among them, ρ0 is the ambient density of the lubricating oil, and D0 represents the coefficient of thermal expansion;

[0065] Calculate the load balance equation for the non - Newtonian fluid thermo - elastohydrodynamic lubrication equation by the separated flow method according to the first characteristic parameter. Among them, the calculation formula for calculating the load balance equation for the non - Newtonian fluid thermo - elastohydrodynamic lubrication equation by the separated flow method according to the first characteristic parameter is:

[0066]

[0067] Among them, F is the total load (N), p h is the oil - film pressure, and p c is the pressure borne by the asperities;

[0068] Calculate the asperity contact pressure according to the Kogut and Etsion contact model of asperity elastic, plastic, and elastoplastic deformations. Among them, the calculation formula for calculating the asperity contact pressure according to the Kogut and Etsion contact model of asperity elastic, plastic, and elastoplastic deformations is:

[0069]

[0070]

[0071] Among them, p c (ω) is the asperity contact pressure, β c is the average radius of the surface rough asperities, K h is the hardness index, v is the Poisson's ratio of the material, H is the hardness of the softer material in the contact surface, is the dimensionless critical deformation amount of the asperities, d * is the dimensionless distance between the average height surface of the asperities and the rigid contact surface, z* is the dimensionless asperity height, E′ is the comprehensive elastic modulus, σ is the root mean square parameter of the surface height, I is the integral operation operator, and α is the exponential value;

[0072] Calculate the normal pressure in the contact area according to the Reynolds equation, the film thickness equation, the viscosity equation, the density equation, the load balance equation, the oil film pressure and the asperity contact pressure. Among them, the calculation formula for calculating the normal pressure in the contact area according to the Reynolds equation, the film thickness equation, the viscosity equation, the density equation, the load balance equation, the oil film pressure and the asperity contact pressure is:

[0073] p a = p h + p c (ω)

[0074] Among them, p a is the normal pressure in the contact area.

[0075] According to some embodiments of the present invention, the numerical calculation of the tooth surface friction force through the non-Newtonian fluid thermo-elastohydrodynamic lubrication equation according to the first characteristic parameter and the true rough surface topography includes:

[0076] Calculate the temperature field control equation of the non-Newtonian fluid thermo-elastohydrodynamic lubrication equation by the separated flow method according to the first characteristic parameter. Among them, the temperature field control equation includes the oil film energy equation, the solid heat conduction equation and the temperature continuity equation. The calculation formula for calculating the temperature field control equation of the non-Newtonian fluid thermo-elastohydrodynamic lubrication equation by the separated flow method according to the first characteristic parameter is:

[0077]

[0078]

[0079]

[0080]

[0081]

[0082] Among them, η* represents the equivalent viscosity, is the dimensionless comprehensive shear stress of the oil film, c f is the specific heat of the lubricating oil, k fis the thermal conductivity of the lubricating oil, u is the flow velocity of the lubricating oil in the x direction, v is the flow velocity of the lubricating oil in the y direction, c1 is the specific heat capacity of gear 1, ρ1 is the density of gear 1, k1 is the thermal conductivity of gear 1, c2 is the specific heat capacity of gear 2, ρ2 is the density of gear 2, k2 is the thermal conductivity of gear 2, u s1 is the flow velocity of the lubricating oil on the surface of gear 1 in the x direction, u s2 is the flow velocity of the lubricating oil on the surface of gear 2 in the x direction, v s1 is the flow velocity of the lubricating oil on the surface of gear 1 in the y direction, v s2 is the flow velocity of the lubricating oil on the surface of gear 2 in the y direction;

[0083] According to the Reynolds equation, the film thickness equation, the viscosity equation, the density equation, the load balance equation and the temperature field control equation, numerical calculations are carried out through dimensionlessization and the Erying rheological model to obtain the oil film shear force. Among them, the calculation formula for obtaining the oil film shear force by numerically calculating according to the Reynolds equation, the film thickness equation, the viscosity equation, the density equation, the load balance equation and the temperature field control equation through dimensionlessization and the Erying rheological model is:

[0084]

[0085]

[0086] Among them, ξ is the dimensionless film thickness direction coordinate, is the dimensionless comprehensive shear stress, G is the shear modulus, is related to the non-Newtonian fluid model, is the dimensionless characteristic shear stress of the lubricating oil, is the dimensionless shear stress of the oil film center layer in the x direction, is the dimensionless shear stress of the oil film center layer in the y direction, τ e is the oil film shear force;

[0087] According to the oil film shear force and the micro-convex body contact pressure, the tooth surface friction force is calculated. Among them, the calculation formula for calculating the tooth surface friction force according to the oil film shear force and the micro-convex body contact pressure is:

[0088] F f =∫∫τ e dxdy+μ c ∫∫p c (ω)dxdy

[0089] Among them, F f is the tooth surface friction force, μ cis the boundary friction coefficient for the preset rough peak friction.

[0090] According to some embodiments of the present invention, the subsurface stress-strain field includes the normal stress in the x-axis direction, the normal stress in the z-axis direction, the normal strain in the x-axis direction, the normal strain in the z-axis direction, the shear stress in the x-z plane, and the shear strain in the x-z plane. Calculating the subsurface stress-strain field according to the normal pressure in the contact area, the residual stress distributed along the tooth surface depth direction, the microhardness, and the tooth surface friction force includes:

[0091] Calculating the subsurface stress-strain field according to the normal pressure in the contact area, the residual stress distributed along the tooth surface depth direction, the microhardness, and the tooth surface friction force by the finite element method.

[0092] According to some embodiments of the present invention, calculating the contact fatigue life of the spiral bevel gear according to the subsurface stress-strain field includes:

[0093] Searching for the critical plane of each material point according to the subsurface stress-strain field by the critical plane method. Wherein, the critical plane of the material point includes the instantaneous normal stress of the material point, the instantaneous normal strain of the material point, and the instantaneous shear strain of the material point. The calculation formula for searching for the critical plane of each material point according to the subsurface stress-strain field by the critical plane method is:

[0094]

[0095] Wherein, σ(t) is the instantaneous normal stress of the material point, ε(t) is the instantaneous normal strain of the material point, γ(t) is the instantaneous shear strain of the material point, θ is the angle between the normal direction of any plane and the positive direction of the x-axis, σ x is the normal stress in the x-axis direction, σ z is the normal stress in the z-axis direction, ε x is the normal strain in the x-axis direction, ε z is the normal strain in the z-axis direction, τ xz is the shear stress in the x-z plane, γ xz is the shear strain in the x-z plane;

[0096] Calculating the maximum shear strain amplitude, the normal strain amplitude, and the mean stress on the critical plane according to the critical plane. Wherein, the calculation formula for calculating the maximum shear strain amplitude, the normal strain amplitude, and the mean stress on the critical plane according to the critical plane is:

[0097]

[0098] Wherein, Δγ max / 2 is the maximum shear strain amplitude, Δε n / 2 is the maximum normal strain amplitude on the critical plane, σ m is the mean stress;

[0099] Calculate the contact fatigue life of the spiral bevel gear based on the Morrow-Brown-Miller fatigue criterion of the critical plane according to the maximum shear strain amplitude, the normal strain amplitude, and the mean stress. Wherein, the calculation formula for calculating the contact fatigue life of the spiral bevel gear based on the Morrow-Brown-Miller fatigue criterion of the critical plane according to the maximum shear strain amplitude, the normal strain amplitude, and the mean stress is:

[0100]

[0101]

[0102] Wherein, E is the elastic modulus, HB is the Brinell hardness, σ′ f is the axial fatigue strength coefficient, ε′ f is the axial fatigue ductility coefficient, b is the fatigue strength index, c is the fatigue ductility index, S is the material constant one, A is the material constant two, B is the material constant three, N f is the contact fatigue life.

[0103] In the second aspect of the present invention, a system for predicting the contact fatigue life of a spiral bevel gear is provided. The system for predicting the contact fatigue life of a spiral bevel gear includes:

[0104] A data acquisition module for acquiring the blank parameters, machining parameters, true rough surface topography, residual stress and microhardness distributed along the tooth surface depth direction of the spiral bevel gear;

[0105] An equivalent meshing characteristic parameter calculation module for calculating equivalent meshing characteristic parameters according to the blank parameters and the machining parameters;

[0106] A first characteristic parameter calculation module for calculating the first characteristic parameter of the bearing contact point on the tooth surface of the bevel gear according to the equivalent meshing characteristic parameters;

[0107] A contact area normal pressure calculation module for numerically calculating the contact area normal pressure through the non-Newtonian fluid thermo-elastohydrodynamic lubrication equation according to the first characteristic parameter and the true rough surface topography;

[0108] A tooth surface friction force calculation module for numerically calculating the tooth surface friction force through the non-Newtonian fluid thermo-elastohydrodynamic lubrication equation according to the first characteristic parameter and the true rough surface topography;

[0109] A subsurface stress-strain field calculation module for calculating the subsurface stress-strain field according to the contact area normal pressure, the residual stress distributed along the tooth surface depth direction, the microhardness, and the tooth surface friction force;

[0110] A contact fatigue life calculation module, configured to calculate the contact fatigue life of the spiral bevel gear according to the subsurface stress-strain field.

[0111] This system obtains the blank parameters, machining parameters, true rough surface topography, residual stress and microhardness distributed along the tooth surface depth direction of the spiral bevel gear; calculates the equivalent meshing characteristic parameters according to the blank parameters and machining parameters; calculates the first characteristic parameters of the load-bearing contact point on the tooth surface of the bevel gear according to the equivalent meshing characteristic parameters; numerically calculates the normal pressure in the contact area through the non-Newtonian fluid thermo-elastohydrodynamic lubrication equation based on the first characteristic parameters and the true rough surface topography; numerically calculates the tooth surface friction force through the non-Newtonian fluid thermo-elastohydrodynamic lubrication equation based on the first characteristic parameters and the true rough surface topography; calculates the subsurface stress-strain field according to the normal pressure in the contact area, the residual stress distributed along the tooth surface depth direction, the microhardness and the tooth surface friction force; calculates the contact fatigue life of the spiral bevel gear according to the subsurface stress-strain field, comprehensively considering the influence of the tooth surface topography-lubrication coupling, residual stress and microhardness on the fatigue performance of the gear, providing a theoretical method and data support for the anti-fatigue design and manufacturing of the spiral bevel gear, and improving the prediction accuracy of the contact fatigue life of the spiral bevel gear.

[0112] In a third aspect of the present invention, there is provided an electronic device for predicting the contact fatigue life of a spiral bevel gear, including at least one control processor and a memory for communicatively connecting with the at least one control processor; the memory stores instructions executable by the at least one control processor, and the instructions are executed by the at least one control processor so that the at least one control processor can execute the above-mentioned method for predicting the contact fatigue life of a spiral bevel gear.

[0113] In a fourth aspect of the present invention, there is provided a computer-readable storage medium storing computer-executable instructions for causing a computer to execute the above-mentioned method for predicting the contact fatigue life of a spiral bevel gear.

[0114] It should be noted that the beneficial effects of the second to fourth aspects of the present invention and the prior art are the same as those of the above-mentioned system for predicting the contact fatigue life of a spiral bevel gear and the prior art, and will not be elaborated here.

[0115] The additional aspects and advantages of the present invention will be partly given in the following description, partly will become obvious from the following description, or will be understood through the practice of the present invention. BRIEF DESCRIPTION OF THE DRAWINGS

[0116] The above and / or additional aspects and advantages of the present invention will become obvious and easy to understand from the description of the embodiments in conjunction with the following drawings, wherein:

[0117] Figure 1 It is a flowchart of a method for predicting the contact fatigue life of spiral bevel gears according to an embodiment of the present invention;

[0118] Figure 2 It is a schematic diagram of a loaded contact analysis model of spiral bevel gears for a method for predicting the contact fatigue life of spiral bevel gears provided by an embodiment of the present invention;

[0119] Figure 3 It is a schematic diagram of the tooth surface contact pressure distribution of the pinion of a spiral bevel gear pair for a method for predicting the contact fatigue life of spiral bevel gears provided by an embodiment of the present invention;

[0120] Figure 4 It is a schematic diagram of the sliding speed, entrainment speed and slide-roll ratio of a spiral bevel gear pair for a method for predicting the contact fatigue life of spiral bevel gears provided by an embodiment of the present invention;

[0121] Figure 5 It is a schematic diagram of the three-dimensional rough tooth surface of a spiral bevel gear measured experimentally for a method for predicting the contact fatigue life of spiral bevel gears provided by an embodiment of the present invention;

[0122] Figure 6 It is a schematic diagram of the residual stress distribution of a spiral bevel gear measured experimentally for a method for predicting the contact fatigue life of spiral bevel gears provided by an embodiment of the present invention;

[0123] Figure 7 It is a schematic diagram of the contact stress and average oil film thickness in the x-axis direction at a certain meshing moment of a spiral bevel gear under two processing methods of grinding and shot peening for a method for predicting the contact fatigue life of spiral bevel gears provided by an embodiment of the present invention;

[0124] Figure 8 It is a schematic diagram of a two-dimensional plane strain elastoplastic finite element model considering lubrication effect and surface integrity for a method for predicting the contact fatigue life of spiral bevel gears provided by an embodiment of the present invention;

[0125] Figure 9 It is a schematic diagram of the time history of each stress component of a material point at a certain depth position during the meshing process of a spiral bevel gear for a method for predicting the contact fatigue life of spiral bevel gears provided by an embodiment of the present invention;

[0126] Figure 10 It is a schematic diagram of the change trend of the contact fatigue life of a spiral bevel gear under two processing methods of grinding and shot peening for a method for predicting the contact fatigue life of spiral bevel gears provided by an embodiment of the present invention;

[0127] Figure 11 It is a schematic diagram of the structure of a system for predicting the contact fatigue life of spiral bevel gears according to an embodiment of the present invention. Detailed implementation manners

[0128] Embodiments of the present invention will be described in detail below. Examples of the embodiments are shown in the accompanying drawings, where like or similar reference numerals denote like or similar elements or elements having like or similar functions throughout. The embodiments described below by referring to the accompanying drawings are exemplary only for explaining the present invention and should not be construed as limiting the present invention.

[0129] In the description of the present invention, if the first, second, etc. are described, it is only for the purpose of distinguishing technical features and should not be construed as indicating or implying relative importance or implicitly indicating the quantity of the indicated technical features or implicitly indicating the sequence of the indicated technical features.

[0130] In the description of the present invention, it should be understood that with regard to the orientation description, such as up, down, etc., the indicated orientation or positional relationship is based on the orientation or positional relationship shown in the accompanying drawings. It is only for the convenience of describing the present invention and simplifying the description, and does not indicate or imply that the indicated device or element must have a specific orientation, be constructed and operated in a specific orientation, and therefore should not be construed as limiting the present invention.

[0131] In the description of the present invention, it should be noted that unless otherwise clearly defined, words such as setting, installation, connection, etc. should be understood in a broad sense, and those skilled in the art can reasonably determine the specific meanings of the above words in the present invention in combination with the specific content of the technical solution.

[0132] Due to its advantages such as strong load-bearing capacity, high contact ratio, and smooth transmission, spiral bevel gears have been widely used in high-speed and heavy-load working conditions in fields such as aerospace, automobiles, and ships. Due to the complex contact geometry, heavy-load working conditions, and tooth surface roughness of the spiral bevel gear pair, the actual lubrication state of the gear is a harsh mixed lubrication state, and the external load is jointly borne by the fluid lubrication area and the micro-convex body contact area. Too low lubricating oil film will cause the direct contact of the tooth surface micro-convex bodies, increasing the stress concentration on the tooth surface; and it is not easy to form a lubricating oil film due to the increase in tooth surface temperature, exacerbating the degree of gear wear. Therefore, studying the mixed lubrication performance of bevel gears has important engineering significance for improving the service life of gears.

[0133] Contact analysis and lubrication analysis are the basis for predicting the contact fatigue life of gears. Due to the complex transmission form and heavy load conditions of spiral bevel gears, the change of the actual contact trajectory during loading has a non-negligible impact on the gear lubrication performance. The prediction of the fatigue life of spiral bevel gears involves the coupling of multiple disciplines such as interfacial contact mechanics, hydrodynamic lubrication mechanics, and fatigue mechanics, and the calculation amount is huge, which brings great challenges to the establishment of a comprehensive and systematic method for predicting the contact fatigue life of spiral bevel gears. On the other hand, traditional gear life assessment methods are based on meeting strength and life theories, lacking a deep understanding of the gear fatigue failure mechanism and the mechanism of the role of surface integrity. Gear surface integrity parameters are crucial for the contact fatigue performance of gears. Therefore, it is urgent to consider the surface integrity and actual load-bearing contact characteristics of spiral bevel gears.

[0134] To solve the above technical defects, referring to Figure 1 , the present invention provides a method for predicting the contact fatigue life of spiral bevel gears, including:

[0135] Step S101, obtaining the blank parameters, machining parameters, true rough surface topography, residual stress and microhardness distributed along the tooth surface depth direction of the spiral bevel gear;

[0136] Step S102, calculating the equivalent meshing characteristic parameters according to the blank parameters and machining parameters;

[0137] Step S103, calculating the first characteristic parameters of the load-bearing contact points on the tooth surface of the bevel gear according to the equivalent meshing characteristic parameters;

[0138] Step S104, performing numerical calculation through the non-Newtonian fluid thermo-elastohydrodynamic lubrication equation according to the first characteristic parameters and the true rough surface topography to obtain the normal pressure in the contact area;

[0139] Step S105, performing numerical calculation through the non-Newtonian fluid thermo-elastohydrodynamic lubrication equation according to the first characteristic parameters and the true rough surface topography to obtain the tooth surface friction force;

[0140] Step S106, calculating the subsurface stress-strain field according to the normal pressure in the contact area, the residual stress distributed along the tooth surface depth direction, the microhardness and the tooth surface friction force;

[0141] Step S107, calculating the contact fatigue life of the spiral bevel gear according to the subsurface stress-strain field.

[0142] This method obtains the blank parameters, machining parameters, true rough surface topography, residual stress and microhardness distributed along the tooth surface depth direction of spiral bevel gears; calculates the equivalent meshing characteristic parameters according to the blank parameters and machining parameters; calculates the first characteristic parameters of the load-bearing contact points on the tooth surface of bevel gears according to the equivalent meshing characteristic parameters; performs numerical calculations through the non-Newtonian fluid thermo-elastohydrodynamic lubrication equation based on the first characteristic parameters and the true rough surface topography to obtain the normal pressure in the contact area; performs numerical calculations through the non-Newtonian fluid thermo-elastohydrodynamic lubrication equation based on the first characteristic parameters and the true rough surface topography to obtain the tooth surface friction force; calculates the subsurface stress-strain field according to the normal pressure in the contact area, the residual stress distributed along the tooth surface depth direction, the microhardness and the tooth surface friction force; calculates the contact fatigue life of spiral bevel gears according to the subsurface stress-strain field, comprehensively considers the influence of tooth surface topography-lubrication coupling, residual stress and microhardness on the fatigue performance of gears, provides a theoretical method and data support for the anti-fatigue design and manufacturing of spiral bevel gears, and improves the prediction accuracy of the contact fatigue life of spiral bevel gears.

[0143] In some embodiments, the equivalent meshing characteristic parameters include the equivalent meshing force, the equivalent meshing moment and the equivalent meshing point position. Calculating the equivalent meshing characteristic parameters according to the blank parameters and machining parameters includes:

[0144] Establish a finite element model according to the blank parameters and machining parameters;

[0145] Perform a loaded contact analysis according to the finite element model to obtain the equivalent meshing force, the equivalent meshing moment and the equivalent meshing point position. Among them, the calculation formula for obtaining the equivalent meshing force according to the finite element model for loaded contact analysis is:

[0146]

[0147] Where, F mf is the equivalent meshing force, N is the number of tooth pairs that are simultaneously engaged at any position of the gear, x is the abscissa of the gear, y is the ordinate of the gear, z is the vertical coordinate of the gear, F mx is the transverse equivalent meshing force, F my is the longitudinal equivalent meshing force, F mz is the vertical equivalent meshing force, F mj，i is the coordinate equivalent meshing force of the i-th tooth pair, F mj is the coordinate equivalent meshing force, and j is the sub-coordinate of the gear;

[0148] The calculation formula for obtaining the equivalent meshing moment according to the finite element model for loaded contact analysis is:

[0149] n j =F mj / F mf (j = x, y, z)

[0150]

[0151] Among them, n j is the direction vector of the equivalent meshing force, and M mj is the equivalent meshing torque, and M mj,i is the equivalent meshing torque of the i-th tooth pair;

[0152] Performing a loaded contact analysis based on the finite element model, the calculation formula for the position of the equivalent meshing point is obtained as follows:

[0153]

[0154] Among them, x m is the abscissa of the position of the equivalent meshing point, y m is the ordinate of the position of the equivalent meshing point, z m is the vertical coordinate of the position of the equivalent meshing point, and the position of the equivalent meshing point is r m (x m , y m , z m ), M mx is the transverse equivalent meshing torque, M my is the longitudinal equivalent meshing torque, and x m，i is the abscissa of the position of the equivalent meshing point of the i-th tooth pair.

[0155] In some embodiments, the first characteristic parameters include the curvature radius of the major axis of the contact ellipse, the curvature radius in the minor axis direction, the relative sliding speed, the entrainment speed and the slide-roll ratio, and the direction angle between the entrainment speed and the minor axis of the contact ellipse. Calculating the first characteristic parameters of the tooth surface load-bearing contact point according to the equivalent meshing characteristic parameters includes:

[0156] Performing curvature calculation according to the position of the equivalent meshing point, the first fundamental form and the second fundamental form of the tooth surface, and obtaining the principal curvatures (k g1 , k g2 , k p1 , k p2 ) of the tooth surfaces of the large and small gears at the load-bearing contact point and the corresponding principal directions (t g1 , t g2 , t p1 , t p2 );

[0157] According to the meshing characteristic parameters, the principal curvatures (k g1 , k g2 , k p1 , k p2 ) and the corresponding principal directions (t g1 , t g2 , t p1 , t p2)Calculate the normal curvature \(k_{1}\) of the large gear in the x-direction at the load-bearing contact point on the tooth surface of the bevel gear, the normal curvature \(k_{2}\) of the small gear, and the induced normal curvature \(k_{I}\). Among them, according to the meshing characteristic parameters, the principal curvatures \((k_{11}, k_{12}, k_{21}, k_{22})\) and the corresponding principal directions \((t_{11}, t_{12}, t_{21}, t_{22})\), the calculation formulas for the normal curvature of the large gear, the normal curvature of the small gear, and the induced normal curvature at the load-bearing contact point on the tooth surface of the bevel gear in the x-direction are as follows: gx , the normal curvature \(k_{2}\) of the small gear px and the induced normal curvature \(k_{I}\) rx , where, according to the meshing characteristic parameters, the principal curvatures \((k_{11}\) g1 , \(k_{12}\) g2 , \(k_{21}\) p1 , \(k_{22}\) p2 ) and the corresponding principal directions \((t_{11}\) g1 , \(t_{12}\) g2 , \(t_{21}\) p1 , \(t_{22}\) p2 ) Calculate the normal curvature of the large gear, the normal curvature of the small gear, and the induced normal curvature at the load-bearing contact point on the tooth surface of the bevel gear in the x-direction. The calculation formula is:

[0158] \(\delta=\arccos(t_{11}\) g1 \(\cdot t_{21}\) p1 / (|t_{11}|\) g1 \(||t_{21}|\)) p1 |))

[0159]

[0160] Among them, is the angle between any tangent line x on the tangent plane at the load-bearing contact point on the tooth surface of the bevel gear and \(t_{11}\) g1 , \(\delta\) is the angle between the \(t_{11}\) direction of the large gear and the \(t_{21}\) direction of the small gear, \(k_{1}\) g1 is the normal curvature of the large gear in the x-direction, \(k_{2}\) p1 is the normal curvature of the small gear, \(k_{I}\) gx is the induced normal curvature; px rx rx Take the maximum and minimum values of the induced normal curvature to obtain the angle \(\theta_{I_{max}}\) corresponding to the maximum value of the induced normal curvature and the angle \(\theta_{I_{min}}\) corresponding to the minimum value of the induced normal curvature. Among them,

[0161] \(\theta_{I_{max}}\) is the angle between the short axis x-axis of the contact ellipse and \(t_{11}\) , \(\theta_{I_{min}}\) is the angle between the long axis y-axis of the contact ellipse and \(t_{11}\) Among them, \(\theta_{I_{max}}\) is the angle between the short axis x-axis of the contact ellipse and \(t_{11}\) g1 , \(\theta_{I_{min}}\) is the angle between the long axis y-axis of the contact ellipse and \(t_{11}\) g1 ;

[0162] Calculate the length of the short semi-axis and the length of the long semi-axis of the contact ellipse according to the equivalent meshing characteristic parameters and the principal curvatures. Among them, the calculation formulas for calculating the length of the short semi-axis and the length of the long semi-axis of the contact ellipse according to the equivalent meshing characteristic parameters and the principal curvatures are:

[0163]

[0164] Among them, a is the length of the short semi-axis of the contact ellipse, b is the length of the short semi-axis of the contact ellipse, \(k_{11}\) aOne of the coefficients of the elliptic integral function, k b Two of the coefficients of the elliptic integral function, F n The normal load, A is a constant one related to the shape of the object, B is a constant two related to the shape of the object, E1 is the elastic modulus of the material of the large gear, E2 is the elastic modulus of the material of the small gear, μ1 is the Poisson's ratio of the material of the large gear, and μ2 is the Poisson's ratio of the material of the small gear;

[0165] Obtain the motion speeds at the load-bearing contact points of the large and small gears according to the basic kinematics of gear transmission. Among them, the calculation formula for obtaining the motion speeds at the load-bearing contact points of the large and small gears according to the basic kinematics of gear transmission is:

[0166] v p = ω p × r p , v g = ω g × r g

[0167] Among them, ω g is the angular velocity vector of the large gear, ω p is the angular velocity vector of the small gear, r g is the radius vector of the large gear, r p is the radius vector of the small gear, v g is the velocity vector of the large gear, v p is the velocity vector of the small gear;

[0168] Calculate the entrainment velocities along the minor and major axes of the contact ellipse according to the velocity vector of the large gear and the velocity vector of the small gear. Among them, the calculation formula for calculating the entrainment velocities along the minor and major axes of the contact ellipse according to the velocity vector of the large gear and the velocity vector of the small gear is:

[0169]

[0170]

[0171] Among them, is the tangential plane component of the velocity vector of the large gear along the minor axis direction of the contact ellipse, is the tangential plane component of the velocity vector of the small gear along the minor axis direction of the contact ellipse, is the component of the velocity vector of the large gear along the tooth surface tangential plane, is the component of the velocity vector of the small gear along the tooth surface tangential plane, is the component of the entrainment velocity along the minor axis direction of the contact ellipse, is the tangential plane component of the entrainment velocity along the major axis direction of the contact ellipse, U e is the entrainment velocity vector, θ e is the angle between the entrainment velocity vector and the minor axis x-axis of the contact ellipse;

[0172] Calculate the entrainment velocity, the angle between the entrainment velocity vector and the minor axis x-axis of the contact ellipse, the sliding velocity, and the angle between the sliding velocity vector and the minor axis x-axis of the contact ellipse according to the entrainment velocities along the minor and major axes of the contact ellipse. The calculation formulas for the entrainment velocity, the angle between the entrainment velocity vector and the minor axis x-axis of the contact ellipse, the sliding velocity, and the angle between the sliding velocity vector and the minor axis x-axis of the contact ellipse are as follows:

[0173]

[0174]

[0175] Wherein, U s is the sliding velocity vector, and θ s is the angle between the sliding velocity vector and the minor axis x-axis of the contact ellipse;

[0176] Calculate the slide-roll ratio according to the entrainment velocity and the sliding velocity. The calculation formula for the slide-roll ratio according to the entrainment velocity and the sliding velocity is as follows:

[0177] AKC = U s / U e

[0178] Wherein, AKC is the slide-roll ratio at the contact point.

[0179] Specifically, by projecting the position of the equivalent meshing point onto the section passing through the gear shaft, the actual load-bearing contact point can be obtained with the help of the tooth surface equation; then, by considering the torsional deformation during gear loading, the relative positions of the two tooth surfaces are adjusted based on the direction vector of the equivalent meshing force so that the contact condition is satisfied at the load-bearing contact point, and then the first characteristic parameter of the load-bearing contact point is solved.

[0180] Specifically, select a hypoid spiral bevel gear for fatigue life prediction. Its main blank parameters and machining parameters are shown in Table 1 and Table 2. Table 1 is the basic design parameters of the spiral bevel gear pair blank, and Table 2 is the gear machining parameters.

[0181] Table 1

[0182]

[0183] Table 2

[0184]

[0185]

[0186] Refer to Figure 2 and Figure 3, based on the blank parameters and machining parameters of spiral bevel gears, a finite element model was established and loaded for contact analysis, and the tooth surface contact pressure distribution of the pinion of the spiral bevel gear pair was solved as shown in Figure 3 shown.

[0187] Referring to Figure 4 , the variation trends of sliding speed, entrainment speed, and slide-roll ratio with the pinion are as shown in Figure 4 shown.

[0188] Referring to Figure 5 and Figure 6 , specifically, two kinds of tooth surfaces of spiral bevel gears processed by grinding and grinding-shot peening were selected for analysis. The three-dimensional topography height information of the tooth surface can be measured by a white light interferometer using the vertical scanning interference method. The tooth surface material can be removed layer by layer by an electrolytic polishing instrument, and the surface residual stress and subsurface residual stress can be measured in combination with an X-ray diffractometer. The microhardness of the gear cross-section can be measured by a Vickers hardness tester.

[0189] In some embodiments, according to the first characteristic parameter and the true rough surface topography, numerical calculations are performed through the non-Newtonian fluid thermo-elastohydrodynamic lubrication equation to obtain the normal pressure in the contact area, including:

[0190] Calculating the Reynolds equation for the non-Newtonian fluid thermo-elastohydrodynamic lubrication equation according to the first characteristic parameter by the separated flow method. Among them, the calculation formula for calculating the Reynolds equation for the non-Newtonian fluid thermo-elastohydrodynamic lubrication equation according to the first characteristic parameter by the separated flow method is:

[0191]

[0192] where p is the pressure, h is the oil film thickness, ρ is the lubricating oil density, η is the lubricating oil viscosity, x is the coordinate in the movement direction, y is the coordinate perpendicular to the movement direction, φ x is the dimensionless flow factor one, φ y is the dimensionless flow factor two;

[0193] Calculating the film thickness equation for the non-Newtonian fluid thermo-elastohydrodynamic lubrication equation according to the first characteristic parameter by the separated flow method. Among them, the calculation formula for calculating the film thickness equation for the non-Newtonian fluid thermo-elastohydrodynamic lubrication equation according to the first characteristic parameter by the separated flow method is:

[0194]

[0195] where h0 is the central film thickness when not deformed, R x is the radius of curvature in the x direction, R y is the radius of curvature in the y direction, E′ is the comprehensive elastic modulus, x′ is the additional coordinate corresponding to x, y′ is the additional coordinate corresponding to y, Ω is the solution domain, and S(x, y) is the roughness height value of the tooth surface;

[0196] Calculate the viscosity equation of the non-Newtonian fluid thermo-elastohydrodynamic lubrication equation by the separated flow method according to the first characteristic parameter. The calculation formula for calculating the viscosity equation of the non-Newtonian fluid thermo-elastohydrodynamic lubrication equation by the separated flow method according to the first characteristic parameter is as follows:

[0197] s0 = β T (T0 - 138) / (lnη0 + 9.67)

[0198]

[0199] Where η0 is the ambient viscosity, z0 is the dimensionless viscosity-pressure index, T is the lubricating oil temperature, T0 is the initial temperature of the lubricating oil, s0 is the dimensionless viscosity-temperature index, and β T is the viscosity-temperature coefficient, and η is the viscosity;

[0200] Calculate the density equation of the non-Newtonian fluid thermo-elastohydrodynamic lubrication equation by the separated flow method according to the first characteristic parameter. The calculation formula for calculating the density equation of the non-Newtonian fluid thermo-elastohydrodynamic lubrication equation by the separated flow method according to the first characteristic parameter is as follows:

[0201]

[0202] Where ρ0 is the ambient density of the lubricating oil, and D0 represents the coefficient of thermal expansion;

[0203] Calculate the load balance equation of the non-Newtonian fluid thermo-elastohydrodynamic lubrication equation by the separated flow method according to the first characteristic parameter. The calculation formula for calculating the load balance equation of the non-Newtonian fluid thermo-elastohydrodynamic lubrication equation by the separated flow method according to the first characteristic parameter is as follows:

[0204]

[0205] Where F is the total load (N), and p h is the oil film pressure, and p c is the pressure borne by the asperities;

[0206] Calculate the asperity contact pressure according to the Kogut and Etsion contact model of asperity elastic, plastic, and elastic-plastic deformation. The calculation formula for calculating the asperity contact pressure according to the Kogut and Etsion contact model of asperity elastic, plastic, and elastic-plastic deformation is as follows:

[0207]

[0208]

[0209] Where p c (ω) is the asperity contact pressure, and β cis the average radius of the surface rough asperities, K h is the hardness index, v is the Poisson's ratio of the material, H is the hardness of the softer material in the contact surface, is the dimensionless critical deformation of the rough peaks, d * is the dimensionless distance between the average height surface of the rough peaks and the rigid contact surface, z * is the dimensionless asperity height, E′ is the comprehensive elastic modulus, σ is the root mean square parameter of the surface height, I is the integral operator, and α is the exponential value;

[0210] The normal pressure in the contact area is calculated according to the Reynolds equation, film thickness equation, viscosity equation, density equation, load balance equation, oil film pressure and micro-convex body contact pressure. Among them, the calculation formula for calculating the normal pressure in the contact area according to the Reynolds equation, film thickness equation, viscosity equation, density equation, load balance equation, oil film pressure and micro-convex body contact pressure is:

[0211] p a = p h + p c (ω)

[0212] Among them, p a is the normal pressure in the contact area.

[0213] In some embodiments, the tooth surface friction force is obtained by numerical calculation according to the first characteristic parameter and the true rough surface topography through the non-Newtonian fluid thermo-elastohydrodynamic lubrication equation, including:

[0214] The temperature field control equation of the non-Newtonian fluid thermo-elastohydrodynamic lubrication equation is calculated according to the first characteristic parameter by the separated flow method. Among them, the temperature field control equation includes the oil film energy equation, solid heat conduction equation and temperature continuity equation. The calculation formula for calculating the temperature field control equation of the non-Newtonian fluid thermo-elastohydrodynamic lubrication equation according to the first characteristic parameter by the separated flow method is:

[0215]

[0216]

[0217]

[0218]

[0219]

[0220] Among them, η* represents the equivalent viscosity, is the dimensionless comprehensive shear stress of the oil film, c f is the specific heat of the lubricating oil, k fis the thermal conductivity of the lubricating oil, u is the flow velocity of the lubricating oil in the x direction, v is the flow velocity of the lubricating oil in the y direction, c1 is the specific heat capacity of gear 1, ρ1 is the density of gear 1, k1 is the thermal conductivity of gear 1, c2 is the specific heat capacity of gear 2, ρ2 is the density of gear 2, k2 is the thermal conductivity of gear 2, u s1 is the flow velocity of the lubricating oil on the surface of gear 1 in the x direction, u s2 is the flow velocity of the lubricating oil on the surface of gear 2 in the x direction, v s1 is the flow velocity of the lubricating oil on the surface of gear 1 in the y direction, v s2 is the flow velocity of the lubricating oil on the surface of gear 2 in the y direction;

[0221] According to the Reynolds equation, film thickness equation, viscosity equation, density equation, load balance equation and temperature field control equation, through dimensionlessization and the Erying rheological model for numerical calculation, the oil film shear force is obtained. Among them, according to the Reynolds equation, film thickness equation, viscosity equation, density equation, load balance equation and temperature field control equation, through dimensionlessization and the Erying rheological model for numerical calculation, the calculation formula for the oil film shear force is:

[0222]

[0223]

[0224] where ξ is the dimensionless film thickness direction coordinate, is the dimensionless comprehensive shear stress, G is the shear modulus, is related to the non-Newtonian fluid model, is the dimensionless characteristic shear stress of the lubricating oil, is the dimensionless shear stress of the oil film center layer in the x direction, is the dimensionless shear stress of the oil film center layer in the y direction, τ e is the oil film shear force;

[0225] According to the oil film shear force and the micro-convex body contact pressure, the tooth surface friction force is calculated. Among them, the calculation formula for calculating the tooth surface friction force according to the oil film shear force and the micro-convex body contact pressure is:

[0226] F f =∫∫τ e dxdy+μ c ∫∫p c (ω)dxdy

[0227] where F f is the tooth surface friction force, μ c is the boundary friction coefficient of the preset rough peak friction.

[0228] Specifically, refer toFigure 7 Taking the rough surface topography of grinding and shot peening as inputs respectively, the pressure distribution and film thickness distribution in the x-axis direction at a certain meshing moment of the non-Newtonian fluid thermo-elastohydrodynamic lubrication of the gear are finally obtained.

[0229] In some embodiments, the subsurface stress-strain field includes the normal stress in the x-axis direction, the normal stress in the z-axis direction, the normal strain in the x-axis direction, the normal strain in the z-axis direction, the shear stress in the x-z plane, and the shear strain in the x-z plane. The subsurface stress-strain field is calculated according to the normal pressure in the contact area, the residual stress distributed along the tooth surface depth direction, the microhardness, and the tooth surface friction force, including:

[0230] The subsurface stress-strain field is calculated by the finite element method according to the normal pressure in the contact area, the residual stress distributed along the tooth surface depth direction, the microhardness, and the tooth surface friction force.

[0231] Specifically, referring to Figure 8 and Figure 9 the time histories of the stress components of the material point at a certain depth position are obtained.

[0232] In some embodiments, the contact fatigue life of the spiral bevel gear is calculated according to the subsurface stress-strain field, including:

[0233] The critical plane of each material point is searched according to the subsurface stress-strain field by the critical plane method. Among them, the critical plane of the material point includes the instantaneous normal stress of the material point, the instantaneous normal strain of the material point, and the instantaneous shear strain of the material point. The calculation formula for searching the critical plane of each material point according to the subsurface stress-strain field by the critical plane method is:

[0234]

[0235] where, σ(t) is the instantaneous normal stress of the material point, ε(t) is the instantaneous normal strain of the material point, γ(t) is the instantaneous shear strain of the material point, θ is the angle between the normal direction of any plane and the positive direction of the x-axis, σ x is the normal stress in the x-axis direction, σ z is the normal stress in the z-axis direction, ε x is the normal strain in the x-axis direction, ε z is the normal strain in the z-axis direction, τ xz is the shear stress in the x-z plane, γ xz is the shear strain in the x-z plane;

[0236] The maximum shear strain amplitude, normal strain amplitude, and mean stress on the critical plane are calculated according to the critical plane. Among them, the calculation formulas for calculating the maximum shear strain amplitude, normal strain amplitude, and mean stress on the critical plane according to the critical plane are:

[0237]

[0238] Among them, Δγ max / 2 is the maximum shear strain amplitude, Δε n / 2 is the maximum normal strain amplitude on the critical plane, and σ m is the mean stress;

[0239] Based on the maximum shear strain amplitude, normal strain amplitude, and mean stress, the contact fatigue life of the spiral bevel gear is calculated according to the Morrow - Brown - Miller fatigue criterion on the critical plane. Among them, the calculation formula for calculating the contact fatigue life of the spiral bevel gear according to the Morrow - Brown - Miller fatigue criterion on the critical plane based on the maximum shear strain amplitude, normal strain amplitude, and mean stress is:

[0240]

[0241]

[0242] Among them, E is the elastic modulus, HB is the Brinell hardness, σ′ f is the axial fatigue strength coefficient, ε′ f is the axial fatigue ductility coefficient, b is the fatigue strength index, c is the fatigue ductility index, S is the material constant one, A is the material constant two, B is the material constant three, and N f is the contact fatigue life.

[0243] Specifically, referring to Figure 10 , the contact fatigue lives of the spiral bevel gears under two processing methods of grinding and shot peening are calculated as shown in Figure 10 .

[0244] Table 3 shows the reference values of other material parameters.

[0245] Table 3

[0246]

[0247] In addition, referring to Figure 11 , an embodiment of the present invention provides a system for predicting the contact fatigue life of a spiral bevel gear, including a data acquisition module 1100, an equivalent meshing feature parameter calculation module 1200, a first feature parameter calculation module 1300, a contact area normal pressure calculation module 1400, a tooth surface friction force calculation module 1500, a subsurface stress - strain field calculation module 1600, and a contact fatigue life calculation module 1700, where:

[0248] The data acquisition module 1100 is used to acquire the blank parameters, processing parameters, true rough surface topography, residual stress and micro - hardness distributed along the tooth surface depth direction of the spiral bevel gear;

[0249] The equivalent meshing characteristic parameter calculation module 1200 is used to calculate equivalent meshing characteristic parameters according to the blank parameters and machining parameters;

[0250] The first characteristic parameter calculation module 1300 is used to calculate the first characteristic parameters of the load-bearing contact points on the tooth surface of the bevel gear according to the meshing characteristic parameters;

[0251] The contact area normal pressure calculation module 1400 is used to perform numerical calculations through the non-Newtonian fluid thermo-elastohydrodynamic lubrication equation according to the first characteristic parameters and the true rough surface topography to obtain the contact area normal pressure;

[0252] The tooth surface friction force calculation module 1500 is used to perform numerical calculations through the non-Newtonian fluid thermo-elastohydrodynamic lubrication equation according to the first characteristic parameters and the true rough surface topography to obtain the tooth surface friction force;

[0253] The subsurface stress-strain field calculation module 1600 is used to calculate the subsurface stress-strain field according to the contact area normal pressure, the residual stress distributed along the tooth surface depth direction, the microhardness and the tooth surface friction force;

[0254] The contact fatigue life calculation module 1700 is used to calculate the contact fatigue life of the spiral bevel gear according to the subsurface stress-strain field.

[0255] This system obtains the blank parameters, machining parameters, true rough surface topography, residual stress and microhardness distributed along the tooth surface depth direction of the spiral bevel gear; calculates equivalent meshing characteristic parameters according to the blank parameters and machining parameters; calculates the first characteristic parameters of the load-bearing contact points on the tooth surface of the bevel gear according to the equivalent meshing characteristic parameters; performs numerical calculations through the non-Newtonian fluid thermo-elastohydrodynamic lubrication equation according to the first characteristic parameters and the true rough surface topography to obtain the contact area normal pressure; performs numerical calculations through the non-Newtonian fluid thermo-elastohydrodynamic lubrication equation according to the first characteristic parameters and the true rough surface topography to obtain the tooth surface friction force; calculates the subsurface stress-strain field according to the contact area normal pressure, the residual stress distributed along the tooth surface depth direction, the microhardness and the tooth surface friction force; calculates the contact fatigue life of the spiral bevel gear according to the subsurface stress-strain field, provides a theoretical method and data support for the anti-fatigue design and manufacture of the spiral bevel gear, and improves the prediction accuracy of the contact fatigue life of the spiral bevel gear.

[0256] It should be noted that the embodiments of this system and the above-mentioned system embodiments are based on the same inventive concept. Therefore, the relevant content of the above method embodiments also applies to the embodiments of this system, and will not be elaborated here.

[0257] The present application also provides an electronic device for predicting the contact fatigue life of spiral bevel gears, including: a memory, a processor, and a computer program stored on the memory and executable on the processor. When the processor executes the computer program, the following is implemented: the method for predicting the contact fatigue life of spiral bevel gears as described above.

[0258] The processor and the memory can be connected through a bus or other means.

[0259] As a non-transitory computer-readable storage medium, the memory can be used to store non-transitory software programs and non-transitory computer-executable programs. In addition, the memory can include high-speed random access memory, and can also include non-transitory memory, such as at least one magnetic disk storage device, a flash memory device, or other non-transitory solid-state storage devices. In some embodiments, the memory can optionally include a memory remotely set relative to the processor, and these remote memories can be connected to the processor through a network. Examples of the above network include but are not limited to the Internet, an enterprise intranet, a local area network, a mobile communication network, and combinations thereof.

[0260] The non-transitory software program and instructions required to implement the method for predicting the contact fatigue life of spiral bevel gears in the above embodiments are stored in the memory. When executed by the processor, the method for predicting the contact fatigue life of spiral bevel gears in the above embodiments is executed. For example, the method steps S101 to S107 described above are executed. Figure 1 in the method steps S101 to S107.

[0261] The present application also provides a computer-readable storage medium storing computer-executable instructions for executing: the method for predicting the contact fatigue life of spiral bevel gears as described above.

[0262] The computer-readable storage medium stores computer-executable instructions, which are executed by a processor or a controller, for example, executed by a processor in the above electronic device embodiment, enabling the above processor to execute the method for predicting the contact fatigue life of spiral bevel gears in the above embodiments. For example, the method steps S101 to S107 described above are executed. Figure 1 in the method steps S101 to S107.

[0263] Those of ordinary skill in the art will understand that all or some of the steps and systems disclosed above can be implemented as software, firmware, hardware, and appropriate combinations thereof. Some or all of the physical components can be implemented as software executed by a processor, such as a central processing unit, a digital signal processor, or a microprocessor, or as hardware, or as an integrated circuit, such as an application specific integrated circuit. Such software can be distributed on a computer-readable medium, which can include a computer storage medium (or non-transitory medium) and a communication medium (or transitory medium). As is well known to those of ordinary skill in the art, the term computer storage medium includes volatile and non-volatile, removable and non-removable media implemented in any method or technology for storing information, such as computer-readable instructions, data structures, program units, or other data. Computer storage media includes, but is not limited to, RAM, ROM, EEPROM, flash memory or other memory technologies, CD-ROM, digital versatile disks (DVD) or other optical disk storage, magnetic cassettes, magnetic tape, magnetic disk storage or other magnetic storage devices, or any other medium that can be used to store the desired information and that can be accessed by a computer. In addition, it is well known to those of ordinary skill in the art that communication media typically contain computer-readable instructions, data structures, program units, or other data in a modulated data signal such as a carrier wave or other transmission mechanism, and can include any information delivery media.

[0264] The embodiments of the present invention have been described in detail above in conjunction with the accompanying drawings, but the present invention is not limited to the above embodiments. Within the scope of knowledge possessed by those of ordinary skill in the art to which the present invention pertains, various changes can be made without departing from the gist of the present invention.

Claims

1. A method for predicting the contact fatigue life of spiral bevel gears, characterized in that, The method for predicting the contact fatigue life of spiral bevel gears includes: Obtaining the blank parameters, machining parameters, true rough surface topography, residual stress and microhardness distributed along the tooth surface depth direction of spiral bevel gears; Calculating equivalent meshing characteristic parameters according to the blank parameters and the machining parameters; Calculating the first characteristic parameters of the load-bearing contact point on the tooth surface of the bevel gear according to the equivalent meshing characteristic parameters; Performing numerical calculation through the non-Newtonian fluid thermo-elastohydrodynamic lubrication equation based on the first characteristic parameters and the true rough surface topography to obtain the microconvex body contact pressure and oil film pressure, and calculating the normal pressure of the contact area according to the microconvex body contact pressure and the oil film pressure. Among them, the calculation formula for calculating the normal pressure of the contact area according to the microconvex body contact pressure and the oil film pressure is: p a = p + p c (ω) where p a is the normal pressure of the contact area, p c (ω) is the asperity contact pressure, and p is the oil film pressure; Performing numerical calculation through the non-Newtonian fluid thermo-elastohydrodynamic lubrication equation based on the first characteristic parameters and the true rough surface topography to obtain the oil film shear force, and calculating the tooth surface friction force according to the oil film shear force and the microconvex body contact pressure. Among them, the calculation formula for calculating the tooth surface friction force according to the oil film shear force and the microconvex body contact pressure is: F f = ∫∫τ e dxdy + μ c ∫∫p c (ω)dxdy Among them, F f is the tooth surface friction force, μ c is the boundary friction coefficient of the preset asperity friction, τ e is the oil film shear force; Calculating the subsurface stress-strain field according to the normal pressure of the contact area, the residual stress distributed along the tooth surface depth direction, the microhardness and the tooth surface friction force; Calculating the contact fatigue life of the spiral bevel gear according to the subsurface stress-strain field.

2. The method for predicting the contact fatigue life of spiral bevel gears according to claim 1, wherein , The equivalent meshing characteristic parameters include equivalent meshing force, equivalent meshing torque and equivalent meshing point position. Calculating the equivalent meshing characteristic parameters according to the blank parameters and the machining parameters includes: Establishing a finite element model according to the blank parameters and the machining parameters; Performing loaded contact analysis according to the finite element model to obtain the equivalent meshing force, the equivalent meshing torque and the equivalent meshing point position. Among them, the calculation formula for obtaining the equivalent meshing force according to the finite element model by performing loaded contact analysis is: Among them, F mf is the equivalent meshing force, N is the number of tooth pairs that are simultaneously meshing at any position of the gear, x is the abscissa of the gear, y is the ordinate of the gear, z is the vertical coordinate of the gear, F mx is the horizontal equivalent meshing force, F my is the longitudinal equivalent meshing force, F mz is the vertical equivalent meshing force, F mj，i is the coordinate equivalent meshing force of the i-th tooth pair, F mj is the coordinate equivalent meshing force, and j is the sub-coordinate of the gear; The calculation formula for obtaining the equivalent meshing torque according to the finite element model by performing loaded contact analysis is: n j = F mj / F mf j = x, y, z where n j is the direction vector of the equivalent meshing force, M mj is the equivalent meshing torque, and M mj,i is the equivalent meshing torque of the i-th tooth pair; The calculation formula for obtaining the equivalent meshing point position according to the finite element model by performing loaded contact analysis is: Among them, x m is the abscissa of the equivalent meshing point position, y m is the ordinate of the equivalent meshing point position, z m is the vertical coordinate of the equivalent meshing point position, and the equivalent meshing point position is r m (x m , y m , z m ), M mx is the transverse equivalent meshing torque, M my is the longitudinal equivalent meshing torque, x m，i is the abscissa of the equivalent meshing point position of the i-th tooth pair.

3. The method for predicting the contact fatigue life of spiral bevel gears according to claim 2, characterized in that, The first characteristic parameters include the curvature radius of the major axis of the contact ellipse, the curvature radius in the minor axis direction, relative sliding speed, entrainment speed and slide-roll ratio, and the direction angle between the entrainment speed and the minor axis of the contact ellipse. Calculating the first characteristic parameters of the load-bearing contact point on the tooth surface of the bevel gear according to the equivalent meshing characteristic parameters includes: Based on the position of the equivalent meshing point, the first fundamental form and the second fundamental form of the tooth surface, curvature calculations are performed to obtain the principal curvatures (k g1 , k g2 , k p1 , k p2 ) of the tooth surfaces of the large and small gears at the load-bearing contact point and the corresponding principal directions (t g1 , t g2 , t p1 , t p2 ); According to the equivalent meshing characteristic parameters, the principal curvatures (k g1 , k g2 , k p1 , k p2 ) and the corresponding principal directions (t g1 , t g2 , t p1 , t p2 ), calculate the normal curvature k gx of the large gear, the normal curvature k px of the small gear and the induced normal curvature k rx of the load-bearing contact point on the tooth surface of the bevel gear in the x-direction. Among them, the calculation formulas for the normal curvature of the large gear, the normal curvature of the small gear and the induced normal curvature of the load-bearing contact point on the tooth surface of the bevel gear in the x-direction according to the equivalent meshing characteristic parameters, the principal curvatures (k g1 , k g2 , k p1 , k p2 ) and the corresponding principal directions (t g1 , t g2 , t p1 , t p2 ) are as follows: δ = arccos(t g1 ·t p1 / (|t g1 ||t p1 |)) Wherein, is the angle between any tangent x and t on the tangent plane of the load-bearing contact point on the tooth surface of the bevel gear, δ is the angle between the t g1 direction of the large gear and the t g1 direction of the small gear, k p1 is the normal curvature of the large gear in the x direction, k gx is the normal curvature of the small gear, k px is the induced normal curvature; rx ​ The induced normal curvature takes maximum and minimum values, and the included angles corresponding to the maximum value of the induced normal curvature and the minimum value of the induced normal curvature are obtained. and the included angles corresponding to the minimum value of the induced normal curvature wherein, is the included angle between the short axis x-axis of the contact ellipse and t g1 ; is the included angle between the long axis y-axis of the contact ellipse and t g1 ; Calculating the minor semi-axis length and major semi-axis length of the contact ellipse according to the equivalent meshing characteristic parameters and the principal curvature. Among them, the calculation formula for calculating the minor semi-axis length and major semi-axis length of the contact ellipse according to the equivalent meshing characteristic parameters and the principal curvature is: where a is the length of the minor semi-axis of the contact ellipse, b is the length of the minor semi-axis of the contact ellipse, k a is the first coefficient of the elliptic integral function, k b is the second coefficient of the elliptic integral function, F n is the normal load, A is the first constant related to the shape of the object, B is the second constant related to the shape of the object, E1 is the elastic modulus of the material of the large gear, E2 is the elastic modulus of the material of the small gear, μ1 is the Poisson's ratio of the material of the large gear, and μ2 is the Poisson's ratio of the material of the small gear; Obtaining the motion speeds at the load-bearing contact points of the large and small gears according to the basic kinematics of gear transmission. Among them, the calculation formula for obtaining the motion speeds at the load-bearing contact points of the large and small gears according to the basic kinematics of gear transmission is: v p = ω p × r p , v g = ω g × r g where ω g is the angular velocity vector of the large gear, ω p is the angular velocity vector of the small gear, r g is the radius vector of the large gear, r p is the radius vector of the small gear, v g is the velocity vector of the large gear, v p is the velocity vector of the small gear; Calculate the entrainment velocity along the minor and major axes of the contact ellipse based on the velocity vectors of the large gear and the small gear. The calculation formula for the entrainment velocity along the minor and major axes of the contact ellipse based on the velocity vectors of the large gear and the small gear is as follows: wherein, is the tangential plane component of the large gear velocity vector along the minor axis direction of the contact ellipse, is the tangential plane component of the small gear velocity vector along the minor axis direction of the contact ellipse, is the component of the large gear velocity vector along the tooth surface tangential plane, is the component of the small gear velocity vector along the tooth surface tangential plane, is the component of the entrainment velocity along the minor axis direction of the contact ellipse, is the tangential plane component of the entrainment velocity along the major axis direction of the contact ellipse, U e is the entrainment velocity vector, θ e is the angle between the entrainment velocity vector and the minor axis x-axis of the contact ellipse; Calculate the entrainment velocity, the angle between the entrainment velocity vector and the x-axis of the minor axis of the contact ellipse, the sliding velocity, and the angle between the sliding velocity vector and the x-axis of the minor axis of the contact ellipse based on the entrainment velocity along the minor and major axes of the contact ellipse. The calculation formula for the entrainment velocity, the angle between the entrainment velocity vector and the x-axis of the minor axis of the contact ellipse, the sliding velocity, and the angle between the sliding velocity vector and the x-axis of the minor axis of the contact ellipse based on the entrainment velocity along the minor and major axes of the contact ellipse is as follows: Among them, U s is the sliding velocity vector, and θ s is the angle between the sliding velocity vector and the minor axis x-axis of the contact ellipse; Calculate the slide-roll ratio based on the entrainment velocity and the sliding velocity. The calculation formula for the slide-roll ratio based on the entrainment velocity and the sliding velocity is as follows: AKC = U s / U e Where AKC is the slide-roll ratio at the contact point.

4. A method for predicting the contact fatigue life of spiral bevel gears according to claim 3, characterized in that Perform numerical calculations on the micro-convex body contact pressure and the oil film pressure through the non-Newtonian fluid thermo-elastohydrodynamic lubrication equation based on the first characteristic parameter and the true rough surface topography, including: Calculate the Reynolds equation for the non-Newtonian fluid thermo-elastohydrodynamic lubrication equation by the separated flow method based on the first characteristic parameter. The calculation formula for calculating the Reynolds equation for the non-Newtonian fluid thermo-elastohydrodynamic lubrication equation by the separated flow method based on the first characteristic parameter is as follows: where p is the oil film pressure, h is the oil film thickness, ρ is the lubricating oil density, η is the lubricating oil viscosity, x is the coordinate in the lateral movement direction, y is the coordinate in the vertical movement direction, and φ x is the dimensionless flow factor one, and φ y is the dimensionless flow factor two; Calculate the film thickness equation for the non-Newtonian fluid thermo-elastohydrodynamic lubrication equation by the separated flow method based on the first characteristic parameter. The calculation formula for calculating the film thickness equation for the non-Newtonian fluid thermo-elastohydrodynamic lubrication equation by the separated flow method based on the first characteristic parameter is as follows: where h0 is the central film thickness when undeformed, R x is the radius of curvature in the x direction, R y is the radius of curvature in the y direction, E′ is the comprehensive elastic modulus, x′ is the additional coordinate corresponding to x, y′ is the additional coordinate corresponding to y, Ω is the solution domain, and S(x, y) is the roughness height value of the tooth surface; Calculate the viscosity equation for the non-Newtonian fluid thermo-elastohydrodynamic lubrication equation by the separated flow method based on the first characteristic parameter. The calculation formula for calculating the viscosity equation for the non-Newtonian fluid thermo-elastohydrodynamic lubrication equation by the separated flow method based on the first characteristic parameter is as follows: s0 = β T (T0 - 138) / (lnη0 + 9.67) Among them, η0 is the environmental viscosity, z0 is the dimensionless viscosity-pressure index, T is the lubricating oil temperature, T0 is the initial temperature of the lubricating oil, s0 is the dimensionless viscosity-temperature index, and β T is the viscosity-temperature coefficient, and η is the viscosity; Calculate the density equation for the non-Newtonian fluid thermo-elastohydrodynamic lubrication equation by the separated flow method based on the first characteristic parameter. The calculation formula for calculating the density equation for the non-Newtonian fluid thermo-elastohydrodynamic lubrication equation by the separated flow method based on the first characteristic parameter is as follows: Where ρ0 is the ambient density of the lubricating oil, and D0 represents the coefficient of thermal expansion; Calculate the load balance equation for the non-Newtonian fluid thermo-elastohydrodynamic lubrication equation by the separated flow method based on the first characteristic parameter. The calculation formula for calculating the load balance equation for the non-Newtonian fluid thermo-elastohydrodynamic lubrication equation by the separated flow method based on the first characteristic parameter is as follows: Among them, F is the total load, p is the oil film pressure, and p c is the pressure borne by the asperity peaks; Calculate the micro-convex body contact pressure according to the Kogut and Etsion contact model for the elastic, plastic, and elastic-plastic deformation of the micro-convex bodies. The calculation formula for calculating the micro-convex body contact pressure according to the Kogut and Etsion contact model for the elastic, plastic, and elastic-plastic deformation of the micro-convex bodies is as follows: Among them, β c is the average radius of the surface rough asperities, K h is the hardness index, H is the hardness of the softer material in the contact surface, is the dimensionless critical deformation amount of the rough peaks, d * is the dimensionless distance between the average height surface of the rough peaks and the rigid contact surface, z * is the dimensionless asperity height, E′ is the comprehensive elastic modulus, σ is the root mean square parameter of the surface height, I is the integral operation operator, and α is the exponential value.

5. A method for predicting the contact fatigue life of spiral bevel gears according to claim 4, characterized in that, Performing numerical calculations on the oil film shear force according to the first characteristic parameter and the true rough surface topography through the non-Newtonian fluid thermo-elastohydrodynamic lubrication equation, including: Calculating the temperature field control equation of the non-Newtonian fluid thermo-elastohydrodynamic lubrication equation by the separated flow method according to the first characteristic parameter, where the temperature field control equation includes the oil film energy equation, the solid heat conduction equation, and the temperature continuity equation. The calculation formula for calculating the temperature field control equation of the non-Newtonian fluid thermo-elastohydrodynamic lubrication equation by the separated flow method according to the first characteristic parameter is: Among them, η* represents the equivalent viscosity, is the dimensionless comprehensive shear stress of the oil film, c f is the specific heat of the lubricating oil, k f is the thermal conductivity of the lubricating oil, u is the flow velocity of the lubricating oil in the x direction, v is the flow velocity of the lubricating oil in the y direction, c1 is the specific heat capacity of gear 1, ρ1 is the density of gear 1, k1 is the thermal conductivity of gear 1, c2 is the specific heat capacity of gear 2, ρ2 is the density of gear 2, k2 is the thermal conductivity of gear 2, u s1 is the flow velocity of the lubricating oil on the surface of gear 1 in the x direction, u s2 is the flow velocity of the lubricating oil on the surface of gear 2 in the x direction, v s1 is the flow velocity of the lubricating oil on the surface of gear 1 in the y direction, v s2 is the flow velocity of the lubricating oil on the surface of gear 2 in the y direction; Performing numerical calculations on the oil film shear force through dimensionlessization and the Eyring rheological model according to the Reynolds equation, the film thickness equation, the viscosity equation, the density equation, the load balance equation, and the temperature field control equation. The calculation formula for performing numerical calculations on the oil film shear force through dimensionlessization and the Eyring rheological model according to the Reynolds equation, the film thickness equation, the viscosity equation, the density equation, the load balance equation, and the temperature field control equation is: where ξ is the dimensionless coordinate in the film thickness direction, is the dimensionless comprehensive shear stress, G is the shear modulus, is related to the non-Newtonian fluid model, is the dimensionless characteristic shear stress of the lubricating oil, is the dimensionless shear stress of the center layer of the oil film in the x direction, is the dimensionless shear stress of the center layer of the oil film in the y direction.

6. A method for predicting the contact fatigue life of spiral bevel gears according to claim 5, characterized in that The subsurface stress and strain field includes the normal stress in the x-axis direction, the normal stress in the z-axis direction, the normal strain in the x-axis direction, the normal strain in the z-axis direction, the shear stress in the x-z plane, and the shear strain in the x-z plane. Calculating the subsurface stress and strain field according to the normal pressure in the contact area, the residual stress distributed along the tooth surface depth direction, the microhardness, and the tooth surface friction force, including: Calculating the subsurface stress and strain field by the finite element method according to the normal pressure in the contact area, the residual stress distributed along the tooth surface depth direction, the microhardness, and the tooth surface friction force.

7. A method for predicting the contact fatigue life of spiral bevel gears according to claim 6, characterized in that, Calculating the contact fatigue life of the spiral bevel gear according to the subsurface stress and strain field, including: Searching for the critical plane of each material point by the critical plane method according to the subsurface stress and strain field, where the critical plane of the material point includes the instantaneous normal stress of the material point, the instantaneous normal strain of the material point, and the instantaneous shear strain of the material point. The calculation formula for searching for the critical plane of each material point by the critical plane method according to the subsurface stress and strain field is: Among them, σ(t) is the instantaneous normal stress of the material point, ε(t) is the instantaneous normal strain of the material point, γ(t) is the instantaneous shear strain of the material point, θ is the angle between the normal direction of any plane and the positive direction of the x-axis, σ x is the normal stress in the x-axis direction, σ z is the normal stress in the z-axis direction, ε x is the normal strain in the x-axis direction, ε z is the normal strain in the z-axis direction, τ xz is the shear stress in the x-z plane, γ xz is the shear strain in the x-z plane; Calculating the maximum shear strain amplitude, the maximum normal strain amplitude, and the mean stress on the critical plane according to the critical plane. The calculation formula for calculating the maximum shear strain amplitude, the maximum normal strain amplitude, and the mean stress on the critical plane according to the critical plane is: where, Δγ max / 2 is the maximum shear strain amplitude, Δε n / 2 is the maximum normal strain amplitude on the critical plane, and σ m is the mean stress; Calculating the contact fatigue life of the spiral bevel gear based on the Morrow-Brown-Miller fatigue criterion of the critical plane according to the maximum shear strain amplitude, the maximum normal strain amplitude, and the mean stress. The calculation formula for calculating the contact fatigue life of the spiral bevel gear based on the Morrow-Brown-Miller fatigue criterion of the critical plane according to the maximum shear strain amplitude, the maximum normal strain amplitude, and the mean stress is: where E is the elastic modulus, HB is the Brinell hardness, σ′ f is the axial fatigue strength coefficient, ε′ f is the axial fatigue ductility coefficient, b is the fatigue strength exponent, c is the fatigue ductility exponent, S is material constant one, A is material constant two, B is material constant three, N f is the contact fatigue life.

8. A spiral bevel gear contact fatigue life prediction system, characterized in that The spiral bevel gear contact fatigue life prediction system includes: A data acquisition module for acquiring blank parameters, machining parameters, true rough surface topography, residual stress and microhardness distributed along the tooth surface depth direction of spiral bevel gears; An equivalent meshing characteristic parameter calculation module for calculating equivalent meshing characteristic parameters according to the blank parameters and the machining parameters; A first characteristic parameter calculation module for calculating the first characteristic parameters of the load-bearing contact points on the tooth surface of bevel gears according to the equivalent meshing characteristic parameters; A contact area normal pressure calculation module for numerically calculating the microconvex body contact pressure and the oil film pressure through the non-Newtonian fluid thermo-elastohydrodynamic lubrication equation based on the first characteristic parameters and the true rough surface topography, and calculating the contact area normal pressure according to the microconvex body contact pressure and the oil film pressure. The calculation formula for calculating the contact area normal pressure according to the microconvex body contact pressure and the oil film pressure is: p a = p + p c (ω) where p a is the normal pressure of the contact area, p c (ω) is the asperity contact pressure, and p is the oil film pressure; A tooth surface friction force calculation module for numerically calculating the oil film shear force through the non-Newtonian fluid thermo-elastohydrodynamic lubrication equation based on the first characteristic parameters and the true rough surface topography, and calculating the tooth surface friction force according to the oil film shear force and the microconvex body contact pressure. The calculation formula for calculating the tooth surface friction force according to the oil film shear force and the microconvex body contact pressure is: F f = ∫∫τ e dxdy + μ c ∫∫p c (ω)dxdy Among them, F f is the tooth surface friction force, μ c is the boundary friction coefficient of the preset asperity friction, τ e is the oil film shear force; A subsurface stress-strain field calculation module for calculating the subsurface stress-strain field according to the contact area normal pressure, the residual stress distributed along the tooth surface depth direction, the microhardness and the tooth surface friction force; A contact fatigue life calculation module for calculating the contact fatigue life of the spiral bevel gear according to the subsurface stress-strain field; 9. A spiral bevel gear contact fatigue life prediction device, characterized in that, Comprising at least one control processor and a memory for communicatingly connecting with the at least one control processor; the memory stores instructions executable by the at least one control processor, and the instructions are executed by the at least one control processor so that the at least one control processor can execute a method for predicting the contact fatigue life of a spiral bevel gear according to any one of claims 1 to 7; 10. A computer-readable storage medium, characterized in that: The computer-readable storage medium stores computer-executable instructions for causing a computer to execute a method for predicting the contact fatigue life of a spiral bevel gear according to any one of claims 1 to 7.

Images