A method for predicting dynamic evolution of cracks in spiral bevel gears

Abstract

The application discloses a kind of spiral bevel gear crack dynamic evolution prediction methods, and the tooth surface distribution force is obtained by bearing tooth contact analysis method, and the crack efficient propagation model of spiral bevel gear is established by finite element method.The stiffness is obtained by meshing simulation to the finite element model containing crack after propagation, and the stiffness is substituted into the dynamics model to solve, and the tooth surface force solved is applied on the tooth surface of finite element model by the method of parameterization programming, to carry out the crack propagation of next stage, and a kind of spiral bevel gear crack dynamic evolution efficient prediction method is established.

Description

TECHNICAL FIELD

[0001] The present application relates to the field of gear system fault dynamics, and particularly relates to a crack dynamic evolution prediction method for spiral bevel gears. BACKGROUND

[0002] Spiral bevel gears are widely used in the fields of aerospace, automobile industry, industrial equipment and the like due to their strong load-carrying capacity, stable transmission and low noise. The crack dynamic evolution research on spiral bevel gears can help to have a deeper understanding of the gear crack failure mechanism, provide a reference for preventing the fatigue crack failure of spiral bevel gear transmission systems, and the crack dynamic evolution of spiral bevel gears is the basis for the design and health assessment of spiral bevel gear transmission systems. Therefore, it is of great significance to research the efficient crack dynamic evolution prediction method for spiral bevel gears.

[0003] A large number of studies have been conducted on the crack propagation of spiral bevel gears in the prior art. In the analysis of the crack propagation of spiral bevel gears, the finite element model is widely used to simulate the crack propagation. Specifically, a geometric model of the spiral bevel gear pair is first established, and then the model is imported into a finite element software. By applying constraints and rotating boundary conditions, a crack is inserted at a specific position to simulate the meshing process of the spiral bevel gear. The direction and size of the crack propagation are determined by the distribution of the stress intensity factor at the crack tip. However, due to the existence of internal and external excitations in the spiral bevel gear transmission system, the meshing force of the spiral bevel gear pair is time-varying, which leads to a large difference between the simulation process of the crack propagation by the finite element method and the real meshing process, and a large error in the simulation results.

[0004] Therefore, there is a need for a new technical solution to solve the above technical problems. SUMMARY

[0005] The technical problem to be solved by the present application is to provide a crack dynamic evolution prediction method for spiral bevel gears, which can be closer to the real simulation of crack propagation characteristics compared with the traditional crack propagation solved by the finite element method, and at the same time improve the calculation efficiency of the fatigue crack propagation.

[0006] To achieve the above-mentioned purpose, the present application can adopt the following technical solutions:

[0007] A crack dynamic evolution prediction method for spiral bevel gears, comprising the following steps:

[0008] (1) The tooth surface points of the pinion and the gear are calculated through the meshing equation and the rotating projection plane equation, and the initial contact points of the pinion and the gear at different meshing moments are calculated through the tooth contact analysis; the potential contact points are calculated through the contact ellipse; the normal load is applied on the tooth surface of the contact surface, the tooth surface deformation is extracted, and the potential contact point flexibility matrix is established; and then the tooth surface node distribution force is solved through the establishment of the deformation coordination equation;

[0009] (2) The geometric model of the spiral bevel gear pair is established by the tooth surface points of the size gears; the distributed force cyclically changing on one tooth surface of the spiral bevel gear is simplified as the distributed force of the tooth surface points at the fixed meshing position; the crack propagation model is established by the geometric model of the spiral bevel gear pair, the simplified distributed force of the tooth surface points at the fixed meshing position and the initial crack inserted at the specified position, the stress intensity factor of the crack tip is calculated through the meshing simulation, the size and direction of the next stage crack propagation are determined, and the shape after the crack propagation is obtained; a single crack growth threshold is set, the crack propagation is carried out, the meshing stiffness of the finite element model of the spiral bevel gear pair containing the crack is calculated through the established load tooth contact analysis model when the crack propagation reaches the single crack growth threshold; the finite element model of the spiral bevel gear pair containing the crack is obtained after the crack propagation simulation based on the crack propagation model;

[0010] (3) The gear transmission system model is discretized and the overall mass matrix M, the overall stiffness matrix K and the overall damping matrix C are calculated; the overall mass matrix M, the overall stiffness matrix K and the overall damping matrix C are integrated in a system to form a dynamic model;

[0011] (4) The meshing stiffness is brought into the dynamic model; the meshing force on the tooth surface is calculated, and then the meshing force is brought into the deformation compatibility equation in the load tooth contact analysis model to obtain new node distributed force; the obtained distributed force is re-substituted into the crack propagation model for simplification, and is simplified as the node distributed force at the fixed position; the size and direction of the next stage crack propagation are calculated through the crack propagation model, so as to form a coupling process until the crack propagation ends.

[0012] Further, in step (1), each initial contact point represents a meshing time and a gear meshing position.

[0013] Further, in step (1), a three-dimensional bilinear isoparametric element is used to establish the finite element model; the unit load normal to the contact point is applied to the tooth surface points of the contact surface of the size gears at different meshing times, and the deformation of the position below a certain distance from the loading point tooth surface is extracted to form the 3D overall flexibility matrix [λ p ] n*i*j and [λ g ] n*i*j , n represents the number of meshing times, i represents the number of tooth surface force application points, and j represents the number of tooth surface displacement extraction points;

[0014] The potential contact point flexibility matrix λ b is:

[0015]

[0016] In the formula, and are the corresponding element values ​​of the flexibility matrix of the potential contact points of the large and small gears, respectively; i represents the contact point from which the displacement is extracted, j represents the contact point from which the force is applied, and n is the number of potential contact points that contact the major axis at this moment of contact;

[0017] And calculate the contact compliance λ c :

[0018] λ c =diag(λ c1 ,λ c2 ,…,λ ci ,…,λ cn ),

[0019] Where E is the elastic modulus; L is the distance between potential contact points; F i is the normal contact force at the ith potential contact point.

[0020] Furthermore, by establishing the deformation coordination equation, the tooth surface node distribution force, that is, the distribution force in the normal direction of the contact point, is solved;

[0021] Among them, the deformation coordination equation is:

[0022]

[0023] Where, F n×1 Represents the distributed force in the normal direction of the contact point; ε n×1 represents the gap of potential contact points; Ste represents the load transfer error; F n_all represents the resultant force of contact forces;

[0024] The normal distribution force of the contact point in a meshing cycle is obtained from the above deformation coordination equation; the point with normal distribution force among the potential contact points is the actual contact point of the spiral bevel gear meshing.

[0025] Furthermore, in step (2), a geometric model of the spiral bevel gear pair is established by three-dimensional software, and the geometric model of the spiral bevel gear pair is imported into finite element software to set boundary conditions, divide the mesh, and extract the node coordinates of the finite element tooth surface;

[0026] The stress of the tooth surface is simplified based on the Paris formula; based on the obtained finite element tooth surface node coordinates and the actual meshing point coordinates obtained by the bearing tooth contact analysis, the closest point between the finite element tooth surface node and the actual meshing point coordinates obtained by the bearing tooth contact analysis is found through programming by using MATLAB software, the potential contact point is matched with the finite element tooth surface node, and the normal distributed force after simplification is applied to the finite element tooth surface node through parameterized programming; the initial crack is introduced into the root of the gear tooth and the meshing simulation is carried out, based on the stress intensity factor theory, the stress intensity factor distribution at the crack tip is calculated, the propagation step is calculated based on the Paris formula, and the propagation direction is calculated based on the MTS theory, so as to carry out crack propagation; the finite element model after crack propagation is substituted into the finite element software to carry out meshing calculation again to obtain the transmission error.

[0027] Further, the meshing stiffness of the finite element model of the crack-containing curved-tooth bevel gear pair is represented as K, and is obtained through the following mathematical function relationship:

[0028]

[0029] STE = θ l × R × sin (δ1) × cos (α t ) × cos (β b )

[0030]

[0031] Wherein, STE is the normal deformation caused by the load, F is the contact force resultant, T is the torque, R is the average pitch, δ1 is the pitch angle, α n is the normal pressure angle, β is the average pitch angle, θ1 is the transmission error angle caused by the elastic deformation, that is, the loaded transmission error angle minus the no-load transmission error angle, α t is the transverse pressure angle, and β b is the helix angle.

[0032] Further, the meshing stiffness K calculated after the crack propagation is substituted into the dynamic model to obtain the required tooth surface meshing force.

[0033] Further, according to the stiffness matrix assembly method in the finite element method, the element node numbers of various types of elements are sequentially assembled into the overall system matrix, and after integrating all the matrices, the overall system dynamic model is obtained:

[0034]

[0035] Wherein, the overall mass matrix M integrates the mass matrix of the spiral bevel gear pair and the shaft coupling mass matrix; the overall stiffness matrix K integrates the meshing stiffness matrix of the spiral bevel gear pair, the shaft coupling stiffness matrix and the bearing support stiffness matrix; the overall damping matrix C is directly calculated using Rayleigh damping; X is the displacement vector, The velocity vector is obtained by taking the first-order derivative of X. The acceleration vector is obtained by taking the second derivative of X.

[0036] Beneficial Effects: This invention obtains the tooth surface distributed force through load-bearing gear tooth contact analysis and establishes an efficient crack propagation model for spiral bevel gears using the finite element method. The expanded cracked finite element model is subjected to meshing simulation to obtain stiffness, which is then substituted into a dynamic model for solution. The solved tooth surface force is then applied to the tooth surface of the finite element model through parametric programming to proceed to the next stage of crack propagation. This establishes an efficient prediction method for the dynamic evolution of cracks in spiral bevel gears.

[0037] Compared with traditional finite element calculations of fatigue crack propagation, the present invention takes into account the characteristics of the mutual coupling between crack propagation and dynamic effects, is closer to the actual fatigue crack propagation process, and the resulting crack propagation changes are more realistic. The crack propagation analyzed based on this model has the advantages of accuracy and reliability, and the calculation efficiency is greatly improved. It can provide a key and reliable foundation for related research such as fault prediction and health assessment of spiral bevel gears.

[0038] The present invention also provides a computer device, comprising a memory, a processor, and a computer program stored in the memory and executable on the processor, wherein the processor implements the steps of the above method when executing the computer program.

[0039] The present invention also provides a computer-readable storage medium having a computer program stored thereon, wherein the computer program implements the steps of the above method when executed by a processor. BRIEF DESCRIPTION OF THE DRAWINGS

[0040] Figure 1 Flowchart of the method for predicting the dynamic evolution of cracks in spiral bevel gears according to the present invention;

[0041] Figure 2 It is a schematic diagram of the gear rotation projection surface;

[0042] Figure 3 Schematic diagram of the gap between corresponding points of conjugate tooth surfaces in the polar coordinate system (LP, θ);

[0043] Figure 4 A schematic diagram of potential contact points;

[0044] Figure 5 Schematic diagram of the distributed force in the normal direction of the contact point during the meshing cycle;

[0045] Figure 6 Schematic diagram of the actual contact points of spiral bevel gear meshing;

[0046] Figure 7 Schematic diagram of the finite element model of the spiral bevel gear;

[0047] Figure 8 Schematic diagram of the normal distribution force applied to the finite element nodes;

[0048] Figure 9 Schematic diagram of meshing plane and transmission error;

[0049] Figure 10 Schematic diagram of the shaft segment beam model;

[0050] Figure 11 Schematic diagram of the bevel gear dynamic model. DETAILED DESCRIPTION

[0051] The present invention will be further described in detail below with reference to the accompanying drawings and specific embodiments.

[0052] Combine Figure 1 As shown, the method for calculating the time-varying meshing stiffness of spiral bevel gears with spalling faults provided by the present invention is applied to the time-varying meshing stiffness calculation of the meshing gears of the spiral bevel gears. The specific steps of the method are as follows.

[0053] 1. Tooth Contact Analysis (TCA)

[0054] (1) Tool rotation surface

[0055] First, obtain the tooth surface point distribution, and obtain the tool rotation surface Σ according to the gear parameters and processing parameters g and Σ p , by vector r g (s g ,θ g ) and r p (s p ,θ p ) can be expressed as:

[0056]

[0057] Where R g Indicates the cutter head radius, s g Indicates the linear blade length, i.e. cutting depth, ɑ g Indicates the tool tooth profile angle, θ g Indicates the angle of tool tip rotation, and the ± symbols correspond to the concave and convex surfaces of the surface obtained by rotating the large gear tool or the small gear tool respectively; R p Indicates the cutter head radius, s p Indicates the linear blade length, i.e. cutting depth, ɑp Indicates the tool tooth profile angle, θ p Indicates the angle of tool tip rotation.

[0058] (2) Coordinates of the tooth surface of the large gear

[0059] Surface Σ g Normal vector n g (θ g ) is expressed as:

[0060]

[0061] Where, the ± signs correspond to the concave normal vector and convex normal vector of the surface obtained by rotating the gear cutter.

[0062] The tooth surface equation (i.e. tooth surface equation) in the wheel blank coordinate system can be obtained by coordinate transformation (Formula (4)):

[0063]

[0064] Where, Indicates the rotation angle of the small wheel blank, M 2g Indicates S g Coordinate transformation matrix to S2. (s2 definition)

[0065] By solving the meshing equation (5) and the rotation projection surface equation (6) simultaneously, the coordinates of any point on the tooth surface of the spiral bevel gear can be obtained.

[0066]

[0067] Where, X m is the x-axis coordinate of the rotation projection surface, Y m is the y-axis coordinate of the rotational projection surface. (r2(z) equation) (parameter definition).

[0068] (3) Pinion tooth surface coordinates

[0069] Similarly, the tool rotation surface is obtained according to the parameters of the pinion and the processing parameters, and then the pinion tooth surface equation is obtained through coordinate transformation. The coordinates of any point on the pinion tooth surface can be obtained by solving the meshing equation and the projection surface equation simultaneously.

[0070] Surface Σ p Normal vector n p (θ p ) is expressed as:

[0071]

[0072] Through coordinate transformation, the tooth surface equation of the pinion is obtained as follows:

[0073]

[0074] The coordinates of any point on the pinion tooth surface can be calculated by simultaneously solving the meshing equation and the rotation projection surface equation:

[0075]

[0076] Where s p Indicates the linear blade length, i.e. cutting depth, θ p Indicates the angle of rotation of the tool tip. It represents the rotation angle of the pinion wheel blank, and the ± symbols correspond to the concave and convex surfaces of the surface obtained by rotating the pinion tool.

[0077] In the gear tooth contact analysis, the initial contact point is calculated, which is the basis for obtaining the actual contact point of the spiral bevel gear meshing through the load-bearing gear tooth contact analysis; in addition, the tooth surface equations of the large and small gears are obtained, and the coordinates of any point on the tooth surface of the large and small gears are calculated, which is the key to obtaining the specific coordinates of the actual contact point.

[0078] 2. Load-bearing Tooth Contact Analysis (LTCA)

[0079] (1) Flexibility matrix and potential contact points

[0080] Based on the gear tooth contact analysis, the load gear tooth contact analysis is carried out. Under load, the meshing state of the spiral bevel gear changes from the contact point to the contact ellipse. A polar coordinate system (L P ,θ), such as Figure 3 As shown, the initial contact point normal vector is moved along the tangent plane by L p , intersecting with the large and small gears at two points respectively. The distance between the two points is the tooth surface clearance d p , calculate the tooth surface clearance under different θ, and find the θ with the minimum clearance min value, this direction is the major axis direction of the contact ellipse, and the maximum gap θ max The value is the direction of the short axis of the contact ellipse, and the point on the long axis of the contact ellipse obtained here is the potential contact point ( Figure 4 ), but note that this is not the actual contact point where the gears mesh.

[0081] TCA can obtain the coordinates of discrete points on the tooth surfaces of large and small gears, and generate the internal mesh of the gears through linear interpolation. A finite element model is established using three-dimensional bilinear isoparametric elements, and constraints other than rotation around the axis are added to the tooth base. At different meshing moments, unit loads in the normal direction of the contact point are applied to the tooth surface points of the large and small gears respectively. The deformation of the tooth surface below 0.2m (m is the module of the gear) at the loading point is extracted and combined into a three-dimensional overall flexibility matrix [λ p ] n*i*j and [λg ] n*i*j , n represents the number of meshing moments, i represents the number of force application points on the tooth surface, and j represents the number of displacement points extracted on the tooth surface.

[0082] According to the overall flexibility matrix of the large and small gears obtained above, the grid points are encrypted by spline interpolation method, and then the flexibility of the interpolation point or grid point closest to the potential contact point is used as the flexibility of the potential contact point to obtain the flexibility matrix λ of the potential contact point p and λ g , and then integrated into the potential contact point flexibility matrix λ according to formula (11) b .

[0083]

[0084] Where, and are the corresponding element values ​​of the flexibility matrix of the potential contact points of the large and small gears respectively; i represents the contact point from which the displacement is extracted, j represents the contact point from which the force is applied, and n is the number of potential contact points that contact the major axis at this contact moment.

[0085] The contact compliance λ is calculated by formula (12): c .

[0086]

[0087] Where E is the elastic modulus; L is the distance between potential contact points; F i is the normal contact force at the ith potential contact point.

[0088] (2) Deformation coordination equation

[0089] Based on the obtained flexibility matrix, deformation coordination iteration is performed according to formula (13).

[0090]

[0091] Where, F n×1 Represents the distributed force in the normal direction of the contact point; ε n×1 represents the gap of potential contact points; Ste represents the load transfer error; F n_all represents the resultant force of the contact force, which is calculated by equation (14).

[0092]

[0093] Where, δ1 represents the cone angle of the spiral bevel gear; ɑ n represents the normal pressure angle of the spiral bevel gear; β represents the average helix angle of the spiral bevel gear.

[0094] (3) Distributed force and actual contact point

[0095] The above deformation coordination equation can be used to obtain the normal distribution force of the contact point in a meshing cycle. The distribution force diagram is shown as follows: Figure 5 As shown, the point with normal distribution force among the potential contact points is the actual contact point of the spiral bevel gear meshing ( Figure 6 ), so the actual contact point position is calculated here through the normal distribution force.

[0096] 3. Establishment of a finite element model for efficient crack propagation

[0097] (1) Based on the coordinates of the discrete points on the tooth surfaces of the large and small gears obtained by the TCA model, a three-dimensional model of the spiral bevel gear is established through three-dimensional software such as UG and Solidworks. The three-dimensional model is imported into finite element software such as ABAQUS and ANSYS to set boundary conditions, divide the mesh, and extract the node coordinates of the finite element tooth surface. The finite element model is as follows: Figure 7 shown.

[0098] (2) The force on the tooth surface is simplified based on the Paris formula. The Paris formula is a crack growth rate model and can be expressed as

[0099] da / dn=C(ΔK) n (15)

[0100] Where a represents the crack extension length, n represents the number of cyclic meshing cycles, C and n are material-related parameters determined by experiments, and the calculation method of ΔK is ΔK = K max -K min , where K max is the maximum value of the stress intensity factor under cyclic loading, where K min is the minimum value of the stress intensity factor under cyclic loading.

[0101] Since the calculation method of ΔK is ΔK=K max -K min , in a complete spiral bevel gear meshing process, K max is unique and K max The corresponding meshing position is the maximum stress position of the tooth root, K min The corresponding position is the moment when the cracked tooth does not participate in meshing and K min =0, so the position with the maximum stress at the tooth root can be used for cyclic loading to replace a complete meshing process calculation, thereby achieving efficient calculation of crack propagation.

[0102] (3) Based on the obtained finite element tooth surface node coordinates and the actual meshing point coordinates obtained by the load-bearing gear tooth contact analysis, the MATLAB software is used to find the point where the finite element tooth surface node is closest to the actual meshing point coordinates obtained by the load-bearing gear tooth contact analysis, and the potential contact point is matched with the finite element tooth surface node. The simplified normal distribution force is applied to the finite element tooth surface node through parametric programming. The schematic diagram of the normal distribution force applied to the finite element node is shown in the figure. Figure 8 shown.

[0103] (4) The initial crack is introduced into the root of the gear tooth and meshing simulation is performed. Based on the stress intensity factor theory, the stress intensity factor distribution at the crack tip is calculated. The expansion step length is calculated based on the Paris formula, and the expansion direction is calculated based on the MTS theory. The crack expansion is then performed. The above process can be achieved using more mature fracture mechanics analysis software currently available on the market, such as FRANC3D.

[0104] Among them, the MTS theory is expressed as:

[0105]

[0106] Where Δθ c is the kink angle, K Ι is the type I stress intensity factor, K ΙΙ is the type II stress intensity factor.

[0107] (5) Calculation of stiffness

[0108] Substituting the finite element model after crack propagation into the finite element software and recalculating the meshing, the transmission error can be obtained, such as Figure 9 As shown, the meshing stiffness of the finite element model can be expressed as K, which can be obtained through the following mathematical function relationship:

[0109]

[0110] STE=θ l ×R×sin(δ1)×cos(α t )×cos(β b ) (18)

[0111]

[0112] where STE is the normal deformation due to the load, F is the resultant contact force, T is the torque, R is the average cone distance, δ1 is the pitch angle, α n is the normal pressure angle, β is the average pitch angle, θ1 is the transmission error angle caused by elastic deformation, that is, the loaded transmission error angle minus the unloaded transmission error angle, α t is the lateral pressure angle, β b is the helix angle.

[0113] 4. Establishment of the kinetic model

[0114] Based on the stiffness calculated after the crack propagation, the stiffness can be substituted into the dynamic model for solution to obtain the required tooth surface meshing force.

[0115] The first step is dynamic modeling: based on the finite element concept, the gear transmission system model is discretized into a series of nodes along the axial direction. These nodes can be divided into shaft segment units, spiral bevel gear pair units and bearing units. By analyzing various units, their differential equations are obtained. After obtaining the mass matrix, stiffness matrix and gyro matrix of the required units, they are integrated into a system to obtain the dynamic equations of the entire gear transmission system and complete the establishment of the dynamic model.

[0116] 1) Shaft segment unit

[0117] Flexible modeling is used to discretize the shaft into nodes along the axial direction. The node locations are generally selected at the ends of the shaft, the cross-section changes, and the gear connections. In order to consider the axial, lateral, and torsional deformations of the shaft, the Timoshenko beam element is used to establish the shaft segment unit dynamic model. The shaft segment unit model is as follows: Figure 10 As shown, a coordinate system O-XYZ is established, where X, Y, and Z are the three displacement directions of nodes A and B, and θ x ,θ y ,θ z are the three rotation directions of nodes A and B. To consider the flexibility of the axis, each axis segment contains two nodes and 12 degrees of freedom:

[0118] u e =[x A y A z A θ xA θ yA θ zA x B y B z B θ xB θ yB θ zB ] T (20)

[0119] 2) Spiral bevel gear pair meshing unit

[0120] A 12-DOF spiral bevel gear meshing unit model is established that takes into account the coupled vibration of various vibration forms such as gear translation, bending, torsion, and pendulum. To simplify the derivation and solution process, the bevel gear is simplified to a rigid cone disk, and the meshing point of the spiral bevel gear pair is set to the midpoint of the tooth width. The dynamic model is as follows Figure 11As shown, local coordinate systems O1-X1Y1Z1 and O2-X2Y2Z2 are established on the driving and driven bevel gears respectively. Figure 11 The axis intersection angle between the driving and driven spiral bevel gears is 90°.

[0121] The lumped mass method is used to establish the dynamic model, simplifying the meshing master and driven spiral bevel gears into two lumped mass points, and using springs to connect the gear pairs to represent the meshing relationship. According to Newton's second law, the dynamic differential equation of the spiral bevel gear meshing unit is established.

[0122] 3) Bearing support unit

[0123] In gear transmission systems, rolling bearings are usually used as support components between the shaft and the housing. To simplify the calculation, the bearing mass and damping are not considered in the modeling, and the time-varying nature of the bearing support stiffness is ignored. The bearing stiffness matrix can be expressed as:

[0124]

[0125] Since most of the coupling directions in the bearing stiffness matrix are relatively small, K b Further simplified to:

[0126] K b =diag{k xx ,k yy ,k zz ,k θxθx ,k θyθy ,0} (22)

[0127] According to the stiffness matrix assembly method in the finite element method, the unit node numbers of various units shown in the above figure are assembled into the overall system matrix in sequence. After integrating all matrices, the dynamic model of the entire system can be obtained:

[0128]

[0129] Wherein, the overall mass matrix M integrates the mass matrix of the spiral bevel gear pair and the shaft coupling mass matrix; the overall stiffness matrix K integrates the meshing stiffness matrix of the spiral bevel gear pair, the shaft coupling stiffness matrix, and the bearing support stiffness matrix; the overall damping matrix C is directly calculated using Rayleigh damping.

[0130] The dynamic force can be calculated by bringing the stiffness obtained by crack extension into the dynamic model.

[0131] 5. Coupling process of efficient crack propagation calculation and dynamic behavior

[0132] First, based on the gear tooth contact analysis model, the initial contact point and the coordinates of the large and small gear tooth surface points are calculated through the gear parameters and processing parameters; based on the load-bearing gear tooth contact analysis model, the potential contact point is calculated, and by applying a normal load on the tooth surface of the contact surface, the tooth surface deformation is extracted and the flexibility matrix is ​​established. Then, by establishing the deformation coordination equation, the actual contact point and the tooth surface node distribution force are solved.

[0133] Secondly, based on the tooth surface point coordinates of the large and small gears obtained from the tooth contact analysis model, a geometric model of the spiral bevel gear was established using 3D software. This geometric model was imported into finite element software, where constraints were added and meshing was performed. After meshing, the tooth surface node coordinates were extracted and matched with the coordinates obtained from the tooth contact analysis model by finding the closest point. Based on the definition of ΔK in the Paris law, the variable tooth surface node distributed force calculated from the loaded tooth contact analysis was simplified to a distributed force for cyclic loading at a fixed position. Based on this simplified distributed force and the tooth surface point coordinates from the finite element, a parametric programming method was used to apply it to the finite element tooth surface, insert an initial crack, and perform crack propagation. This established an efficient crack propagation model. A single crack growth threshold was set. When the crack reached this single crack growth threshold, the stiffness was calculated based on the gear model after crack propagation.

[0134] Finally, based on the obtained stiffness of the spiral bevel gear and the established rigid-flexible coupling dynamic model, the stiffness is incorporated into the dynamic model for calculation to obtain the vibration response of the gear transmission system. The solved dynamic force is substituted into the deformation coordination equation in the load-bearing gear tooth contact analysis model to calculate the distributed force on the gear tooth surface for a new meshing cycle. The new distributed force is introduced into the established efficient crack propagation model and further simplified into a fixed-position distributed force. The simplified force is then applied to the tooth surface of the finite element model that has undergone crack propagation, and a new crack propagation calculation is continued, forming a coupling process between efficient crack propagation calculation and dynamic behavior until the crack propagation ends. In this way, an efficient prediction method for the dynamic evolution of cracks in spiral bevel gears is established.

[0135] The present invention has many specific application paths, and the above description is only a preferred embodiment of the present invention. It should be noted that those skilled in the art can make several improvements and modifications without departing from the principles of the present invention, and such improvements and modifications should also be considered as the scope of protection of the present invention.

Claims

1. A method for predicting the dynamic evolution of cracks in spiral bevel gears, characterized in that: The following steps are involved: (1) The tooth surface points of the large and small gears are calculated by the meshing equation and the rotational projection surface equation, and the initial contact points of the large and small gears at different meshing moments are calculated by tooth contact analysis; Calculate potential contact points using contact ellipses; By applying a normal load on the contact surface, the tooth surface deformation is extracted and the potential contact point flexibility matrix is ​​established; then, by establishing the deformation coordination equation, the tooth surface node distribution force is solved; (2) A geometric model of a bevel gear pair is established using the tooth surface points of the large and small gears; the cyclically changing distributed force on one tooth surface of the bevel gear is simplified into a tooth surface point distributed force at a fixed meshing position; a crack propagation model is established using the geometric model of the bevel gear pair, the simplified tooth surface point distributed force at a fixed meshing position, and an initial crack inserted at a specified position; and the stress intensity factor at the crack tip is calculated through meshing simulation; Determine the size and direction of the crack expansion in the next stage and obtain the shape of the crack after expansion; A single crack growth threshold is set and crack propagation is performed. When the crack propagates to the single crack growth threshold, the meshing stiffness of a finite element model of a cracked spiral bevel gear pair is calculated using the established load-bearing gear tooth contact analysis model. The cracked spiral bevel gear pair finite element model is a finite element model obtained after crack propagation simulation based on the crack propagation model. A geometric model of the spiral bevel gear pair is established using three-dimensional software, and the geometric model is imported into the finite element software to set boundary conditions, divide the mesh, and extract the node coordinates of the finite element tooth surface. The force on the tooth surface is simplified based on the Paris formula. Based on the obtained finite element tooth surface node coordinates and the actual meshing point coordinates obtained from the load-bearing gear tooth contact analysis, MATLAB software is used to program to find the point closest to the finite element tooth surface node and the actual meshing point coordinates obtained from the load-bearing gear tooth contact analysis. The potential contact point is matched with the finite element tooth surface node, and the simplified normal distribution force is applied to the finite element tooth surface node through parametric programming. The initial crack is introduced into the root of the gear tooth and a meshing simulation is performed. Based on the stress intensity factor theory, the stress intensity factor distribution at the crack tip is calculated. The expansion step size is calculated based on the Paris formula, and the expansion direction is calculated based on the MTS theory, thereby performing crack expansion. Substitute the finite element model after crack propagation into the finite element software to recalculate the meshing to obtain the transmission error; The meshing stiffness of the finite element model of the cracked spiral bevel gear pair is expressed as K, which is obtained by the following mathematical function relationship: STE=θ l ×R×sin(δ1)×cos(α t )×cos(β b ) where STE is the normal deformation due to the load, F is the resultant contact force, T is the torque, R is the average cone distance, δ1 is the pitch angle, α n is the normal pressure angle, β is the average pitch angle, θ1 is the transmission error angle caused by elastic deformation, that is, the loaded transmission error angle minus the unloaded transmission error angle, α t is the lateral pressure angle, β b is the helix angle; (3) Discretizing the gear transmission system model and calculating the overall mass matrix M, the overall stiffness matrix K, and the overall damping matrix C; integrating the overall mass matrix M, the overall stiffness matrix K, and the overall damping matrix C into a system to form a dynamic model; (4) The meshing stiffness is brought into the dynamic model; the tooth surface meshing force is calculated, and then the meshing force is brought into the deformation coordination equation in the load-bearing gear tooth contact analysis model to obtain a new node distribution force, and the obtained distribution force is substituted back into the crack propagation model for simplification, simplified to a node distribution force at a fixed position, and the crack propagation size and direction of the next stage are calculated by the crack propagation model, thereby forming a coupling process until the crack propagation ends.

2. The method for predicting the dynamic evolution of cracks in spiral bevel gears according to claim 1, wherein: In step (1), each initial contact point represents a meshing moment and a gear meshing position.

3. The method for predicting the dynamic evolution of cracks in spiral bevel gears according to claim 1, wherein: In step (1), a finite element model is established using a three-dimensional bilinear isoparametric unit; at different meshing moments, a unit load in the normal direction of the contact point is applied to the tooth surface points of the large and small gears respectively, and the deformation at a certain distance below the loading point is extracted to form a three-dimensional overall flexibility matrix [λ p ] n*i*j and [λ g ] n*i*j , n represents the number of meshing moments, i represents the number of force application points on the tooth surface, and j represents the number of displacement points extracted on the tooth surface; Potential contact point flexibility matrix λ b for: Where, and are the corresponding element values ​​of the flexibility matrix of the potential contact points of the large and small gears, respectively; i represents the contact point from which the displacement is extracted, j represents the contact point from which the force is applied, and n is the number of potential contact points that contact the major axis at this moment of contact; And calculate the contact compliance λ c : l c =diag(λ c1 ,l c2 ,…,l ci ,…,l cn ), Where E is the elastic modulus; L is the distance between potential contact points; F i is the normal contact force at the ith potential contact point.

4. The method for predicting the dynamic evolution of cracks in spiral bevel gears according to claim 1, wherein: Then, by establishing the deformation coordination equation, the tooth surface node distribution force, that is, the distribution force in the normal direction of the contact point, is solved; Among them, the deformation coordination equation is: Where, F n×1 Represents the distributed force in the normal direction of the contact point; ε n×1 represents the gap of potential contact points; Ste represents the load transfer error; F n_all represents the resultant force of contact forces; The normal distribution force of the contact point in a meshing cycle is obtained from the above deformation coordination equation; the point with normal distribution force among the potential contact points is the actual contact point of the spiral bevel gear meshing.

5. The method for predicting the dynamic evolution of cracks in spiral bevel gears according to claim 1, wherein: The meshing stiffness K calculated based on the above crack expansion is substituted into the dynamic model to solve and obtain the required tooth surface meshing force.

6. The method for predicting the dynamic evolution of cracks in spiral bevel gears according to claim 5, characterized in that: According to the stiffness matrix assembly method in the finite element method, the unit node numbers of various units are assembled into the overall system matrix in sequence. After integrating all matrices, the dynamic model of the entire system is obtained: Wherein, the overall mass matrix M integrates the mass matrix of the spiral bevel gear pair and the shaft system coupling mass matrix; The overall stiffness matrix K integrates the meshing stiffness matrix of the spiral bevel gear pair, the shaft system coupling stiffness matrix and the bearing support stiffness matrix; The overall damping matrix C is directly calculated using Rayleigh damping; X is the displacement vector, The velocity vector is obtained by taking the first-order derivative of X. The acceleration vector is obtained by taking the second derivative of X.

7. A computer device comprising a memory, a processor, and a computer program stored in the memory and executable on the processor, wherein: When the processor executes the computer program, the steps of the method according to any one of claims 1 to 6 are implemented.

8. A computer-readable storage medium having a computer program stored thereon, characterized in that: When the computer program is executed by a processor, the steps of the method according to any one of claims 1 to 6 are implemented.

Images