Contact stiffness calculation and dynamic modeling method for non-Hertzian contact of planetary gear train

Abstract

The present invention discloses a contact stiffness calculation and dynamic modeling method for non-Hertzian contact of a planetary gear train. The steps of the contact stiffness calculation method are as follows: obtaining gear tooth surface test data; using an involute weighted compensation method to match the gear tooth surface test data with the gear tooth surface curvature to calculate the contact stiffness. The dynamic modeling method calculates the contact stiffness according to the above method; then the axial compressive stiffness, bending stiffness, shear stiffness, calculated contact stiffness, etc. are brought into the planetary gear train dynamics model, and the ode45 calculation method is used to solve the dynamic equation to obtain the gear tooth vibration data; the gear tooth vibration data is analyzed to obtain the force on the gear tooth. The contact form of the gear tooth surface is considered as non-Hertzian contact, the microscopic morphology of the gear tooth surface is fully considered, the calculation error of the contact stiffness is reduced, the accuracy of the dynamic modeling is improved, and the accuracy of the contact characteristic analysis, vibration characteristic estimation, power loss calculation, and life prediction is improved.

Description

Technical Field

[0001] The invention relates to the technical field of simulation analysis, and in particular to a contact stiffness calculation and dynamic modeling method for non-Hertzian contact of a planetary gear train. Background Art

[0002] The planetary gear train has a large transmission ratio, strong load-bearing capacity, and a variety of transmission forms. The transmission ratios of different transmission forms are different. For example, under the same gear tooth parameters and sun gear input speed conditions, the transmission ratios of the planetary carrier output and the ring gear output will be different if the transmission form is different. The planetary gear train has a variety of transmission forms, so it is widely used in high-speed, heavy-load or multi-speed transmission mechanisms. Dynamic modeling of the planetary gear train can analyze its motion characteristics, calculate power loss, etc., and can be used to predict the life of the planetary gear train.

[0003] The current dynamic modeling of planetary gear trains does not fully consider the microscopic morphology of the gear tooth surface, and considers the contact form of the gear tooth surface as Hertzian contact, which causes large calculation errors and reduces the accuracy of further contact characteristic analysis, vibration characteristic estimation, power loss calculation and even life prediction. Summary of the invention

[0004] The purpose of the present invention is to overcome the deficiency in the prior art that the above-mentioned contact form is considered as Hertzian contact and the micromorphology of the gear tooth surface is not fully considered, and to provide a contact stiffness calculation and dynamic modeling method for non-Hertzian contact of a planetary gear system. By fully considering the micromorphology of the gear tooth surface and considering the contact form of the gear tooth surface as non-Hertzian contact, the calculation error is reduced and the accuracy of subsequent calculation simulation is improved.

[0005] In order to achieve the above-mentioned object of the invention, the present invention provides the following technical solutions:

[0006] In a first aspect, an embodiment of the present invention provides a method for calculating contact stiffness of non-Hertzian contact of a planetary gear train, comprising the following steps:

[0007] S1, obtaining gear tooth surface test data; using the involute weighted compensation method to match the gear tooth surface test data with the gear tooth surface curvature, and calculating the contact stiffness.

[0008] Preferably, the step S1 specifically includes the following steps:

[0009] S11, microscopically scanning the gear tooth surface by using an ultra-depth-of-field microscopic scanner to obtain gear tooth surface test data;

[0010] S12, compensating the gear tooth surface test data and calculating the fractal dimension of the planetary gear tooth surface profile;

[0011] The involute weighted processing method is adopted to combine the pressure angle of different contact positions on the gear tooth surface and the angle between the photographic direction of the microscanner and the center of the scanned gear tooth to obtain the compensation ratio of different contact positions; the gear tooth surface test data is compensated, and the compensated gear tooth surface test data is recorded as the first gear tooth surface test data;

[0012] The least squares midline slope is calculated after taking the double logarithm of the power spectrum function; the fractal dimension is obtained according to the least squares midline slope;

[0013] S13, select points for the gear surface test data, and perform curve fitting on the selected points to obtain the contour layer data fitting curve:

[0014] The section corresponding to the maximum value of the gear surface test data is called the surface layer, and the contour line of the surface layer is recorded as the characteristic meshing line; on the characteristic meshing line, with the data peak of the characteristic meshing line as the center, points are taken at equal intervals on the left and right to select points for gear surface test data; curve fitting is performed on the selected points to obtain the contour layer data fitting curve;

[0015] S14, calculating the actual contact area and the contact area of ​​the fitting curve by fitting the profile layer data curve;

[0016] S15, calculating the surface contact coefficient, and performing segmented calculation of the dimensionless normal load according to the deformation stage type and the value range of the fractal dimension;

[0017] S16, calculate the non-Hertzian contact stiffness.

[0018] Preferably, the step S12 calculates the fractal dimension D according to the following formula:

[0019]

[0020] Where k is the least squares midline slope calculated by taking the double logarithm of the power spectrum function;

[0021] The step S12 compensates the gear surface test data according to the formula shown below:

[0022] K=[tan(a i )-a i ]'

[0023] ε i =K / tan(θ)

[0024] H i =ε i *H

[0025] Among them, α i is the pressure angle corresponding to the meshing position, K is the slope of the involute function corresponding to the pressure angle, ε iis the compensation ratio of the contact position, θ is the angle between the photographic direction of the microscanner and the center of the scanned gear tooth, H is the height deviation of the microscan of the gear tooth surface corresponding to the meshing position, H i is the height after compensation.

[0026] In the step S13, the points are taken at equal intervals on the left and right sides, and the distance n between the points is calculated according to the following formula:

[0027] n=N / D

[0028] Where N is the total number of points of the characteristic curve and D is the fractal dimension.

[0029] Preferably, the step S14 calculates the actual contact area according to the following formula:

[0030] According to the critical contact point area a c With the maximum contact area a l The actual contact area A is calculated by segmentation r , calculated as follows:

[0031] when a l >a c hour:

[0032]

[0033] when a l c hour:

[0034]

[0035] Where D is the fractal dimension;

[0036] According to the elastic-plastic critical deformation boundary condition, the critical contact area a between the plastic deformation and elastic deformation of a single micro-convex body is c The calculation formula is as follows:

[0037]

[0038] Where G is the roughness amplitude of the original contact rough surface; b is related to the hardness I and yield strength σ of the softer material. y The relevant coefficient is calculated as q is the yield strength of the softer material σ y The coefficient related to the comprehensive elastic modulus E is calculated as

[0039] Among them, the maximum contact area a l ​According to the non-Hertzian contact fractal theory, the distribution function of the contact size of the micro-convex body with a contact area of ​​a is calculated; the distribution function is as follows:

[0040]

[0041] Among them, n is the distance between the points, D is the fractal dimension, λ c is the original surface contact coefficient, which is selected according to the gear surface accuracy requirements, a l is the area of ​​the maximum contact point between the bonding surfaces, recorded as the maximum contact area.

[0042] Preferably, the surface contact coefficient λ in step S15 c The calculation formula is as follows:

[0043]

[0044] Where F is the external load, ΣS represents the sum of the surface areas of the two contact surfaces, ΣS = 2πB1(R1±R2) (“+” represents external contact, “-” represents internal contact, R1 and R2 are the curvature radii of cylinder 1 and cylinder 2 respectively, and B1 is the effective tooth width); x h represents the comprehensive curvature radius coefficient, B is the effective contact length of the contact object, here B = 20; F is the normal load between the surfaces; ω(B, E) represents the coefficient related to the effective contact length B and the comprehensive elastic modulus, C1, C2 represent coefficients related to the surface geometry, where C1 = 1, C2 = 0.5 when the cylindrical surfaces are in contact;

[0045] Step S15 performs segmented calculation of the dimensionless normal load according to the deformation stage type and the value range of the fractal dimension D. The calculation method is as follows:

[0046] Critical contact area a c With the maximum contact area a l The relationship is a c l When , the mechanical interface contact is in the elastic-plastic deformation stage, and the dimensionless normal load is calculated as follows:

[0047] When D≠1.5, the dimensionless normal load F * for:

[0048]

[0049] When D = 1.5, the dimensionless normal load F * for:

[0050] ​

[0051] In the formula, K φ is the contact deformation stiffness of the gear tooth surface, g1(D), g2(D) are the simplified calculation formulas for the fractal dimension, g1(D), g2(D), and The calculation is as follows:

[0052]

[0053] a c >a l When the mechanical interface contact is in the complete plastic deformation stage, the dimensionless normal load is calculated as follows:

[0054]

[0055] Preferably, the step S16 calculates the normal non-Hertzian contact stiffness between the mechanical joint surfaces according to the following formula:

[0056]

[0057] in, is the dimensionless normal contact stiffness; E is the comprehensive elastic modulus of the contact object surface, Wherein, μ1 and μ2 refer to the Poisson's ratio of contact objects 1 and 2, respectively, E1 and E2 refer to the elastic modulus of contact objects 1 and 2, respectively; D is the fractal dimension of the bonding surface; a l is the area of ​​the maximum contact point between the bonding surfaces, a c is the critical area between plastic deformation and elastic deformation of a single micro-convex body; is the dimensionless true contact area between the contact surfaces.

[0058] In a second aspect, an embodiment of the present invention provides a method for dynamic modeling of non-Hertzian contact of a planetary gear train, comprising the following steps:

[0059] S1, calculating the contact stiffness of the non-Hertzian contact according to any of the above-mentioned contact stiffness calculation methods for the non-Hertzian contact of the planetary gear train;

[0060] S2, calculate the axial compressive stiffness, bending stiffness and shear stiffness according to the gear meshing energy method;

[0061] S3, bring the contact stiffness, axial compressive stiffness, bending stiffness and shear stiffness into the translational torsional dynamics model of the planetary gear train, use the ode45 calculation method to solve the dynamic equations and obtain the vibration data of the gear teeth;

[0062] S4, analyzing the vibration data of the gear teeth, and obtaining the force of the gear teeth by analyzing the force and deformation of the gear teeth.

[0063] Preferably, the step S3 calculates the comprehensive meshing stiffness k according to the following formula: e :

[0064]

[0065] where k b is the bending stiffness, k s is the shear stiffness, k a is the axial compressive stiffness, k f is the matrix stiffness, is the non-Hertzian contact stiffness, where 1 and 2 are the numbers of the meshing gears.

[0066] Preferably, the planetary gear train dynamics equation in step S3 is expressed in the form of a matrix as shown below:

[0067] Mq″+[c b +c e ]q′+[k b +k e ]q=T(t)+F(t)

[0068] Where M is the mass matrix, k b / c b is the support stiffness matrix, k e / c e is the meshing stiffness matrix, T is the external excitation force, and F is the internal excitation force.

[0069] Compared with the prior art, the present invention has the following beneficial effects:

[0070] 1. Consider the contact form of the gear tooth surface as non-Hertzian contact, fully consider the microscopic morphology of the gear tooth surface, reduce the calculation error of the contact stiffness, improve the accuracy of dynamic modeling, and improve the accuracy of contact characteristic analysis, vibration characteristic estimation, power loss calculation, and life prediction.

[0071] 2. The ultra-depth-of-field microscopic test data of the gear surface is brought into the gear contact stiffness calculation process, replacing the Hertz contact stiffness of the traditional planetary gear dynamics modeling, and the calculation error is reduced.

[0072] 3. The involute weighted compensation method is used to restore the original curvature of the gear tooth surface to solve the matching problem between the curvature of the gear tooth surface and the microscopic scanning.

[0073] 4. Dimensional point sampling to collect surface data can collect surface data more comprehensively to obtain surface data that better represents the characteristics of the gear teeth. BRIEF DESCRIPTION OF THE DRAWINGS

[0074] Figure 1 Flow chart of a method for calculating contact stiffness of non-Hertzian contact of a planetary gear train according to exemplary embodiment 1 of the present invention;

[0075] Figure 2 It is a schematic diagram of a profile layer data fitting curve of exemplary embodiment 1 of the present invention;

[0076] Figure 3 This is a flow chart of a method for dynamic modeling of non-Hertzian contact of a planetary gear train according to exemplary embodiment 2 of the present invention. DETAILED DESCRIPTION

[0077] The present invention is further described in detail below in conjunction with test examples and specific implementation methods. However, this should not be understood as the scope of the above subject matter of the present invention being limited to the following embodiments, and all technologies realized based on the content of the present invention belong to the scope of the present invention.

[0078] Example 1

[0079] This embodiment provides a method for calculating contact stiffness of non-Hertzian contact of a planetary gear train, comprising the following steps:

[0080] S1, obtaining gear tooth surface test data; using the involute weighted compensation method to match the gear tooth surface test data with the gear tooth surface curvature, and calculating the contact stiffness.

[0081] In this embodiment, a three-dimensional ultra-depth microscanner is used to obtain gear tooth surface test data; the gear tooth surface test data obtained by the ultra-depth microscanner is brought into the gear tooth contact stiffness calculation process, replacing the Hertz contact stiffness of the traditional planetary gear dynamics modeling, and the involute weighted compensation method is used to solve the matching problem between the gear tooth surface curvature and the microscanning data. The contact stiffness calculation method described in this embodiment fully considers the microscopic morphology of the gear tooth surface, considers the contact form of the gear tooth surface as non-Hertz contact, reduces the calculation error, and improves the accuracy of subsequent calculation simulation.

[0082] Among them, Figure 1 As shown, step S1 specifically includes the following steps:

[0083] S11, microscopically scanning the gear tooth surface by using an ultra-depth-of-field microscopic scanner to obtain gear tooth surface test data;

[0084] S12, compensating the gear tooth surface test data and calculating the fractal dimension of the planetary gear tooth surface profile;

[0085] The involute weighted processing method is adopted to combine the pressure angle of different contact positions on the gear tooth surface and the angle between the photographic direction of the microscanner and the center of the scanned gear tooth to obtain the compensation ratio of different contact positions; the gear tooth surface test data is compensated, and the compensated gear tooth surface test data is recorded as the first gear tooth surface test data;

[0086] Taking double logarithm of the power spectrum function of the gear tooth surface test data, and then calculating the least squares midline slope, the least squares midline slope is used to calculate the fractal dimension;

[0087] S13, select points for the gear surface test data, and perform curve fitting on the selected points to obtain the contour layer data fitting curve:

[0088] The section corresponding to the maximum value of the original measured gear surface test data is called the surface layer, and the contour line of the surface layer is recorded as the characteristic meshing line; on the characteristic meshing line, with the data peak of the characteristic meshing line as the center, points are taken at equal intervals on the left and right to select points for gear surface test data; curve fitting is performed on the selected points to obtain the contour layer data fitting curve;

[0089] S14, calculating the actual contact area and the contact area of ​​the fitting curve by fitting the profile layer data curve;

[0090] S15, calculating the surface contact coefficient, and performing segmented calculation of the dimensionless normal load according to the deformation stage type and the value range of the fractal dimension;

[0091] S16, calculate the non-Hertzian contact stiffness.

[0092] The contact stiffness calculation method described in this embodiment fully considers the microscopic morphology of the gear tooth surface and considers the contact form of the gear tooth surface as non-Hertzian contact, thereby reducing calculation errors and improving the accuracy of subsequent calculation simulations.

[0093] Exemplarily, step S12 calculates the fractal dimension D according to the formula shown in formula (1):

[0094]

[0095] Wherein, k is the least squares midline slope calculated by taking the double logarithm of the power spectrum function of the test waveform signal on the gear tooth surface;

[0096] The rough surface refers to the surface formed by the gear surface test data; the surface profile curve refers to the curve formed by taking points on the gear surface data on the characteristic meshing line and fitting it, which is the profile layer data fitting curve described later. The fractal dimension of the rough surface and the surface profile curve is calculated according to formula (1) respectively, and the fractal dimension of the rough surface D is obtained. s and the fractal dimension D of the surface contour curve p , the relationship between the two is shown in formula (2):

[0097] D s =D p +1 (2)

[0098] In the actual calculation process, the fractal dimension D of the rough surface is taken for the rough surface. s , take the fractal dimension D of the surface contour curve for the surface contour curve p In the subsequent calculation process of this embodiment, since the corresponding calculation is performed for the rough surface, the fractal dimension D used in the following calculation process is the fractal dimension D of the rough surface. s .

[0099] Exemplarily, step S12 compensates the gear tooth surface test data according to formulas (3) to (5):

[0100] K=[tan(a i )-a i ]' (3)

[0101] ε i =K / tan(θ) (4)

[0102] H i =ε i *H (5)

[0103] Among them, α i is the pressure angle corresponding to the meshing position, K is the slope of the involute function corresponding to the pressure angle, ε i is the compensation ratio of the contact position, θ is the angle between the photographic direction of the microscanner and the center of the scanned gear tooth, H is the height deviation of the microscan of the gear tooth surface corresponding to the meshing position, H i is the height after compensation.

[0104] Exemplarily, when selecting points in step S13, the point distance of the selected points is calculated according to the formula shown in formula (6):

[0105] The section corresponding to the maximum value of the original measured gear surface test data is called the surface layer, and the contour line of the surface layer is recorded as the characteristic meshing line; on the characteristic meshing line, with the data peak of the characteristic meshing line as the center, points are taken at equal intervals on the left and right, and the point number distance n is:

[0106] n=N / D (6)

[0107] Where N is the total number of points of the characteristic meshing line, and D is the fractal dimension.

[0108] After selecting the points, all the selected points form a set A, and the selected points are fitted with a curve by the quadratic fitting method to obtain the contour layer data fitting curve, such as Figure 2 shown.

[0109] Exemplarily, step S14 calculates the actual contact area according to the formulas shown in equations (7) to (10):

[0110] According to the elastic-plastic critical deformation boundary condition, the critical contact area a between the plastic deformation and elastic deformation of a single micro-convex body is c The calculation formula is:

[0111]

[0112] Where G is the roughness amplitude of the original contact rough surface; b is the hardness I and yield strength σ of the softer material. y The relevant coefficient is calculated as q is the yield strength of the softer material σ y The coefficient related to the comprehensive elastic modulus E is calculated as

[0113] According to the non-Hertzian contact fractal theory, the contact size of a micro-convex body with a contact area of ​​a obeys the distribution function:

[0114]

[0115] Among them, n is the distance between the points, D is the fractal dimension, λ c is the surface contact coefficient, which is selected according to the gear surface accuracy requirements, a l is the area of ​​the maximum contact point between the bonding surfaces, recorded as the maximum contact area;

[0116] The contour layer data fitting curve is as follows Figure 2 As shown, the maximum contact area a l That is Figure 2 The shaded part of the fitting curve in the figure; the maximum contact area a is calculated according to formula (8) l ;

[0117] According to the critical contact point area a c Calculate the actual contact area A r :

[0118] when a l >a c hour:

[0119]

[0120] when a l c hour:

[0121]

[0122] ​Exemplarily, the step S15 calculates the surface contact coefficient according to the formula shown in formula (11):

[0123] According to the surface contact characteristics, the surface contact coefficient is constructed under the condition of satisfying the Hertz contact assumption. The surface contact coefficient λ c The calculation formula is:

[0124]

[0125] In the formula, F e is the external load, i.e. the normal load between the surfaces; ∑S represents the sum of the surface areas of the two contact surfaces, ∑S = 2πB1(R1±R2) (“+” represents external contact, “-” represents internal contact, R1 and R2 are the curvature radii of cylinder 1 and cylinder 2 respectively, and B1 is the effective tooth width); x h represents the comprehensive curvature radius coefficient, B is the effective contact length of the contact object, and here B=20; ω(B,E) represents the coefficient related to the effective contact length B and the comprehensive elastic modulus, C1 and C2 represent coefficients related to the surface geometry. Here, when the cylindrical surfaces are in contact, the values ​​are C1=1 and C2=0.5.

[0126] Exemplarily, step S15 calculates the dimensionless normal load according to the formulas shown in equations (12) to (15):

[0127] The dimensionless normal load between the bonding surfaces can be classified and calculated according to the value of the fractal dimension D. The formula is:

[0128] 1) Elastic-plastic deformation stage, there is a c l ;

[0129] When D≠1.5, the dimensionless normal load is

[0130]

[0131] When D = 1.5, the dimensionless normal load is

[0132]

[0133] In the formula, K φ is the contact deformation stiffness of the gear tooth surface, which is selected as K according to the gear material. φ =5.731×10 8 N / m, A r is the actual contact area, A a is the theoretical contact area, is the actual area contact coefficient, ​is the contact coefficient of the critical contact point area, where g1(D), g2(D), The calculation method of is as follows:

[0134]

[0135] 2) When the mechanical interface contact is in the complete plastic deformation stage, there is a c >a l , the dimensionless normal load is calculated as:

[0136]

[0137] Exemplarily, step S16 calculates the normal non-Hertzian contact stiffness between the mechanical joint surfaces according to the formula shown in formula (16):

[0138]

[0139] in, is the normal non-Hertzian contact stiffness; E is the comprehensive elastic modulus of the contact object surface, Wherein, μ1 and μ2 refer to the Poisson's ratio of contact objects 1 and 2, respectively, E1 and E2 refer to the elastic modulus of contact objects 1 and 2, respectively; D is the fractal dimension of the bonding surface; a l is the area of ​​the maximum contact point between the bonding surfaces, a c is the critical contact area between the plastic deformation and elastic deformation of a single micro-convex body; is the true contact area between dimensionless contact surfaces, that is, the actual area contact coefficient.

[0140] Example 2

[0141] like Figure 3 As shown, this embodiment provides a method for dynamic modeling of non-Hertzian contact of a planetary gear train, comprising the following steps:

[0142] S1, calculating the contact stiffness of non-Hertzian contact according to the method described in Example 1;

[0143] S2, calculate the axial compressive stiffness, bending stiffness and shear stiffness according to the gear meshing energy method;

[0144] S3, bring the contact stiffness, axial compressive stiffness, bending stiffness and shear stiffness into the translational torsional dynamics model of the planetary gear train, use the ode45 calculation method to solve the dynamic equations and obtain the vibration data of the gear teeth;

[0145] S4, analyzing the vibration data of the gear teeth, and obtaining the force of the gear teeth by analyzing the force and deformation of the gear teeth.

[0146] Exemplarily, step S3 calculates the comprehensive meshing stiffness k according to the formula shown in formula (17): e :

[0147]

[0148] where k b is the bending stiffness, k s is the shear stiffness, k a is the axial compressive stiffness, k f is the matrix stiffness, is the non-Hertzian contact stiffness, where 1 and 2 are the numbers of the meshing gears.

[0149] Exemplarily, the planetary gear train dynamics equation in step S3 is expressed in the form of a matrix as shown in formula (18):

[0150] Mq″+[c b +c e ]q′+[k be +k e ]q=T(t)+F(t)(18)

[0151] Where M is the mass matrix, q represents the yield strength of the softer material, and c b is the support damping matrix, k be / c b is the support stiffness matrix, c e is the meshing damping matrix, k e / c e is the meshing stiffness matrix, T is the external excitation force, and F is the internal excitation force.

[0152] The above is only a detailed description of the specific implementation of the present invention, rather than a limitation of the present invention. Various substitutions, modifications and improvements made by those skilled in the relevant art without departing from the principle and scope of the present invention should be included in the protection scope of the present invention.

Claims

1. A method for calculating contact stiffness of non-Hertzian contact of a planetary gear train, characterized in that: The following steps are involved: S1, obtaining gear tooth surface test data; using the involute weighted compensation method to match the gear tooth surface test data with the gear tooth surface curvature, and calculating the contact stiffness; The step S1 specifically includes the following steps: S11, microscopically scanning the gear tooth surface by using an ultra-depth-of-field microscopic scanner to obtain gear tooth surface test data; S12, compensating the gear tooth surface test data and calculating the fractal dimension of the planetary gear tooth surface profile; The involute weighted processing method is adopted to combine the pressure angle of different contact positions on the gear tooth surface and the angle between the photographic direction of the microscanner and the center of the scanned gear tooth to obtain the compensation ratio of different contact positions; the gear tooth surface test data is compensated, and the compensated gear tooth surface test data is recorded as the first gear tooth surface test data; The least squares midline slope is calculated after taking the double logarithm of the power spectrum function; the fractal dimension is obtained according to the least squares midline slope; S13, select points for the gear surface test data, and perform curve fitting on the selected points to obtain the contour layer data fitting curve: The section corresponding to the maximum value of the gear surface test data is called the surface layer, and the contour line of the surface layer is recorded as the characteristic meshing line; on the characteristic meshing line, with the data peak of the characteristic meshing line as the center, points are taken at equal intervals on the left and right to select points for gear surface test data; curve fitting is performed on the selected points to obtain the contour layer data fitting curve; S14, calculating the actual contact area and the contact area of ​​the fitting curve by fitting the profile layer data curve; S15, calculating the surface contact coefficient, and performing segmented calculation of the dimensionless normal load according to the deformation stage type and the value range of the fractal dimension; S16, calculate the non-Hertzian contact stiffness.

2. The method for calculating contact stiffness of non-Hertzian contact of a planetary gear train according to claim 1, characterized in that: The step S12 calculates the fractal dimension D according to the formula shown below: Where k is the least squares midline slope calculated by taking the double logarithm of the power spectrum function; The step S12 compensates the gear surface test data according to the formula shown below: K=[tan(a i )-to i ]′ e i =K / tan(θ) H i =ε i *H Among them, α i is the pressure angle corresponding to the meshing position, K is the slope of the involute function corresponding to the pressure angle, ε i is the compensation ratio of the contact position, θ is the angle between the photographic direction of the microscanner and the center of the scanned gear tooth, H is the height deviation of the microscan of the gear tooth surface corresponding to the meshing position, H i is the height after compensation.

3. The method for calculating contact stiffness of non-Hertzian contact of a planetary gear train according to claim 1, characterized in that: In step S13, points are taken at equal intervals on the left and right sides, and the point number distance n of the taken points is calculated according to the formula shown below: n=N / D Where N is the total number of points of the characteristic curve and D is the fractal dimension.

4. The method for calculating contact stiffness of non-Hertzian contact of a planetary gear train according to claim 1, characterized in that: The step S14 calculates the actual contact area according to the following formula: According to the critical contact point area a c With the maximum contact area a l The actual contact area A is calculated by segmentation r , calculated as follows: when a l >a c hour: when a l c hour:​ Where D is the fractal dimension; According to the elastic-plastic critical deformation boundary condition, the critical contact area a between the plastic deformation and elastic deformation of a single micro-convex body is c The calculation formula is as follows: Where G is the roughness amplitude of the original contact rough surface; b is related to the hardness I and yield strength σ of the softer material. y The relevant coefficient is calculated as q is the yield strength of the softer material σ y The coefficient related to the comprehensive elastic modulus E is calculated as Among them, the maximum contact area a l According to the non-Hertzian contact fractal theory, the distribution function of the contact size of the micro-convex body with a contact area of ​​a is calculated; the distribution function is as follows: Among them, n is the distance between the points, D is the fractal dimension, λ c is the original surface contact coefficient, which is selected according to the gear surface accuracy requirements, a l is the area of ​​the maximum contact point between the bonding surfaces, recorded as the maximum contact area.

5. The method for calculating contact stiffness of non-Hertzian contact of a planetary gear train according to claim 1, characterized in that: The surface contact coefficient λ in step S15 c The calculation formula is as follows: Where F is the external load, that is, the normal load between the surfaces, ∑S represents the sum of the surface areas of the two contact surfaces, ∑S = 2πB1(R1±R2), "+" represents external contact, "-" represents internal contact, R1 and R2 are the curvature radii of cylinder 1 and cylinder 2 respectively, and B1 is the effective tooth width; x h represents the comprehensive curvature radius coefficient, B is the effective contact length of the contact object, and here B=20; ω(B,E) represents the coefficient related to the effective contact length B and the comprehensive elastic modulus, C1, C2 represent coefficients related to the geometric shape of the curved surface. Here, when the cylindrical surface is in contact, the values ​​are C1 = 1, C2 = 0.5; E represents the comprehensive elastic modulus; Step S15 performs segmented calculation of the dimensionless normal load according to the deformation stage type and the value range of the fractal dimension D. The calculation method is as follows: Critical contact area a c With the maximum contact area a l The relationship is a c l When , the mechanical interface contact is in the elastic-plastic deformation stage, and the dimensionless normal load is calculated as follows:​ When D≠1.5, the dimensionless normal load F * for: When D = 1.5, the dimensionless normal load F * for: In the formula, K φ is the contact deformation stiffness of the gear tooth surface, g1(D), g2(D) are the simplified calculation formulas for the fractal dimension, g1(D), g2(D), and The calculation is as follows: a c >a l When the mechanical interface contact is in the complete plastic deformation stage, the dimensionless normal load is calculated as follows: A r Indicates the actual contact area; A α represents the theoretical contact area, represents the actual area contact coefficient, is the contact coefficient of the critical contact point area.

6. The method for calculating contact stiffness of non-Hertzian contact of a planetary gear train according to claim 1, characterized in that: The step S16 calculates the normal non-Hertzian contact stiffness between the mechanical joint surfaces according to the formula shown below: in, is the dimensionless normal contact stiffness; E is the comprehensive elastic modulus of the contact object surface, Wherein, μ1 and μ2 refer to the Poisson's ratio of contact objects 1 and 2, respectively, E1 and E2 refer to the elastic modulus of contact objects 1 and 2, respectively; D is the fractal dimension of the bonding surface; a l is the area of ​​the maximum contact point between the bonding surfaces, a c is the critical area between plastic deformation and elastic deformation of a single micro-convex body; is the real contact area between dimensionless contact surfaces; A r Represents the actual contact area; λ c Represents the surface contact coefficient.

7. A method for dynamic modeling of non-Hertzian contact of a planetary gear train, characterized in that: The following steps are involved: S1, calculating the contact stiffness of non-Hertzian contact according to the method described in any one of claims 1 to 6; S2, calculate the axial compressive stiffness, bending stiffness and shear stiffness according to the gear meshing energy method; S3, bring the contact stiffness, axial compressive stiffness, bending stiffness and shear stiffness into the translational torsional dynamics model of the planetary gear train, use the ode45 calculation method to solve the dynamic equations and obtain the vibration data of the gear teeth; S4, analyzing the vibration data of the gear teeth, and obtaining the force of the gear teeth by analyzing the force and deformation of the gear teeth.

8. The method for dynamic modeling of non-Hertzian contact of a planetary gear train according to claim 7, characterized in that: The step S3 calculates the comprehensive meshing stiffness k according to the formula shown below: e : where k b is the bending stiffness, k s is the shear stiffness, k a is the axial compressive stiffness, k f is the matrix stiffness, is the non-Hertzian contact stiffness, where 1 and 2 are the numbers of the meshing gears.

9. The method for dynamic modeling of non-Hertzian contact of a planetary gear train according to claim 7, characterized in that: The planetary gear train dynamics equation in step S3 is expressed in the form of a matrix as shown below: Where M is the mass matrix, k b / c b is the support stiffness matrix, k e / c e is the meshing stiffness matrix, T is the external excitation force, and F is the internal excitation force; c b represents the support damping matrix, c e represents the mesh damping matrix.

Images