Spiral bevel gear crack dynamic evolution prediction method

Abstract

The invention discloses a spiral bevel gear crack dynamic evolution prediction method, and the method comprises the steps: obtaining a tooth surface distribution force through a bearing gear tooth contact analysis method, and building a spiral bevel gear crack efficient propagation model through a finite element method. Carrying out meshing simulation on the expanded finite element model containing the cracks to obtain rigidity, substituting the rigidity into the dynamic model for solving, applying the solved tooth surface force to the tooth surface of the finite element model through a parameterized programming method, and carrying out crack expansion of the next stage, and a spiral bevel gear crack dynamic evolution efficient prediction method is established.

Description

Technical Field

[0001] The invention relates to the field of gear system fault dynamics, and in particular to a method for predicting the dynamic evolution of cracks in a spiral bevel gear. Background Art

[0002] Spiral bevel gears are widely used in aerospace, automotive, industrial equipment and other fields due to their strong load-bearing capacity, smooth transmission and low noise. The study of the dynamic evolution of cracks in spiral bevel gears can provide a deeper understanding of the mechanism of gear crack failure and provide a reference for preventing fatigue crack failures in spiral bevel gear transmission systems. The dynamic evolution of cracks in spiral bevel gears is the basis for the design and health assessment of spiral bevel gear transmission systems, so it is of great significance to study the efficient prediction method of the dynamic evolution of cracks in spiral bevel gears.

[0003] A lot of research has been done on the crack propagation of spiral bevel gears in the prior art. When analyzing the crack propagation of spiral bevel gears, the widely used method is to establish a finite element model for crack propagation. Specifically, the geometric model of the spiral bevel gear pair is first established, and then the model is imported into the finite element software. By applying constraints and rotation boundary conditions, cracks are inserted at specific positions to simulate the meshing process of the spiral bevel gear, and the direction and size of crack propagation are determined by the distribution of the stress intensity factor at the crack tip. However, during the meshing process, 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, resulting in a large difference between the process of simulating crack propagation using the finite element method and the actual meshing process, and the simulation results will produce large errors.

[0004] Therefore, a new technical solution is needed to solve the above technical problems. Summary of the invention

[0005] The technical problem to be solved by the present invention is to provide a method for predicting the dynamic evolution of cracks in spiral bevel gears, which can simulate crack propagation characteristics closer to reality than the traditional crack propagation solved by finite element method, and at the same time improve the calculation efficiency of fatigue crack propagation.

[0006] To achieve the above object, the present invention can adopt the following technical solutions:

[0007] A method for predicting the dynamic evolution of cracks in a spiral bevel gear comprises the following steps:

[0008] (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 the gear tooth contact analysis; the potential contact points are calculated by the contact ellipse; the tooth surface deformation is extracted and the potential contact point flexibility matrix is ​​established by applying a normal load on the tooth surface of the contact surface; and the tooth surface node distribution force is solved by establishing the deformation coordination equation;

[0009] (2) Establishing a geometric model of an arc bevel gear pair through the tooth surface points of large and small gears; simplifying the cyclically changing distributed force on one tooth surface of the arc bevel gear into the tooth surface point distributed force at a fixed meshing position; establishing a crack propagation model through the geometric model of the arc bevel gear pair, the simplified tooth surface point distributed force at a fixed meshing position and inserting an initial crack at a specified position, and calculating the stress intensity factor of the crack tip through meshing simulation; determining the size and direction of the crack propagation in the next stage, and obtaining the shape after crack propagation; setting a single crack growth threshold, and performing crack propagation. When the crack propagates to the single crack growth threshold, the meshing stiffness of the finite element model of the arc bevel gear pair containing cracks is calculated through the established load-bearing gear tooth contact analysis model; the finite element model of the arc bevel gear pair containing cracks is a finite element model obtained after crack propagation simulation based on the crack propagation model;

[0010] (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;

[0011] (4) The meshing stiffness is introduced into the dynamic model; the meshing force of the tooth surface is calculated, and then the meshing force is introduced 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 re-substituted 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 through the crack propagation model, thereby forming a coupling process until the crack propagation ends.

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

[0013] Furthermore, in step (1), a finite element model is established using a three-dimensional bilinear isoparametric unit; 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, and the deformation at a position below a certain distance from the tooth surface of 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;

[0014] Potential contact point flexibility matrix λ b for:

[0015]

[0016] In the formula, and are the corresponding element values ​​of the compliance 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;

[0017] And calculate the contact flexibility λ c :

[0018]

[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] In the formula, 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 meshing of the spiral bevel gear.

[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 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 by the load-bearing gear tooth contact analysis, the MATLAB software is used to program 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 initial crack is introduced into the root of the gear tooth and the meshing simulation is performed. Based on the stress intensity factor theory, the stress intensity factor distribution at the crack tip is calculated, the extension step length is calculated based on the Paris formula, and the extension direction is calculated based on the MTS theory, so as to perform crack extension; the finite element model after crack extension is substituted into the finite element software for re-meshing calculation to obtain the transmission error.

[0027] Furthermore, the meshing stiffness of the finite element model of the spiral bevel gear pair with cracks is expressed as K, which is obtained by the following mathematical function relationship:

[0028]

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

[0030]

[0031] 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.

[0032] Furthermore, the meshing stiffness K calculated based on the above crack extension is substituted into the dynamic model to solve and obtain the required tooth surface meshing force.

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

[0034]

[0035] In the formula, 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 second derivative of X is the acceleration vector.

[0036] Beneficial effect: The present invention obtains the tooth surface distribution force through the load-bearing gear tooth contact analysis method, and establishes an efficient crack expansion model for arc bevel gears through the finite element method. The expanded finite element model containing cracks is subjected to meshing simulation to obtain the stiffness, which is substituted into the dynamic model for solution. The solved tooth surface force is applied to the tooth surface of the finite element model through the parametric programming method to carry out the next stage of crack expansion, and an efficient prediction method for the dynamic evolution of cracks in arc bevel gears is established.

[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 on which a computer program is stored, and the computer program implements the steps of the above method when executed by a processor. BRIEF DESCRIPTION OF THE DRAWINGS

[0040] Figure 1 It is a flow chart of the method for predicting the dynamic evolution of cracks in spiral bevel gears of the present invention;

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

[0042] Figure 3 It is a 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 point of spiral bevel gear meshing;

[0046] Figure 7 It is a 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] Fig. 9 Schematic diagram of meshing plane and transmission error;

[0049] Fig.10 It is a schematic diagram of the shaft segment beam model;

[0050] Fig.11 Schematic diagram of bevel gear dynamics model. DETAILED DESCRIPTION

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

[0052] Combination 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 calculation of the time-varying meshing stiffness of the large and small gears of the spiral bevel gears meshing with each other. The specific steps of the method are as follows.

[0053] 1. Tooth Contact Analysis (TCA)

[0054] (1) Tool rotation surface

[0055] First, the tooth surface point distribution is obtained, and the surface Σ of the tool rotation is obtained 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] In the formula, R g Indicates the radius of the cutter head, s g Indicates the linear blade length, i.e. the cutting depth, ɑ g Indicates the tool tooth profile angle, θ g Indicates the angle of tool tip rotation, and the ± signs correspond to the concave and convex surfaces of the surface obtained by rotating the large gear tool or the small gear tool; R p Indicates the radius of the cutter head, s p Indicates the linear blade length, i.e. the 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 The 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 large gear cutter.

[0062] Through coordinate transformation (Formula (4)), the tooth surface envelope equation (i.e. tooth surface equation) in the wheel blank coordinate system can be obtained: r2(s g ,θg , ).

[0063]

[0064] In the formula, 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 rotation projection surface. (r2(z) equation) (parameter definition).

[0068] (3) Pinion tooth surface coordinates

[0069] Similarly, the tool rotation surface is obtained according to the pinion parameters and 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 The 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 obtained by combining the meshing equation and the rotational projection surface equation:

[0075]

[0076] In the formula, s p Indicates the linear blade length, i.e. the cutting depth, θ p Indicates the angle of rotation of the tool tip. It represents the rotation angle of the pinion wheel blank, and the ± signs 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, so as to calculate the coordinates of any point on the tooth surface of the large and small gears, 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, and 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 This direction is the major axis direction of the contact ellipse, and the maximum gap θ max The value is the direction of the minor axis of the contact ellipse. The point on the major 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 surface of large and small gears, and generate the internal mesh of the gears through linear interpolation. The finite element model is established using three-dimensional bilinear isoparametric units, and constraints other than rotation around the axis are added to the tooth base; unit loads in the normal direction of the contact point are applied to the tooth surface points of the large and small gears at different meshing moments, and the deformation of the tooth surface below 0.2m (m is the module of the gear) at the loading point is extracted, and the 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 gears obtained above, the grid points are encrypted by spline interpolation, 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] 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 where the displacement is extracted, j represents the contact point where 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 flexibility λ 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] In the formula, 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 schematic diagram of the distribution force is shown in Figure 5 As shown in Figure 2, 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 by the normal distributed force.

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

[0097] (1) Based on the coordinates of the discrete points on the gear tooth surfaces 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 position of maximum stress at 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 extension step length is calculated based on the Paris formula. The extension direction is calculated based on the MTS theory, and the crack extension is 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] In the formula, Δθ 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 Fig. 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 extension, the stiffness can be substituted into the dynamic model for solution to obtain the required tooth surface meshing force.

[0115] The first 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, arc bevel gear pair units and bearing units. The differential equations of various units are obtained by analyzing them. 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 positions 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: Fig.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 pair meshing unit model is established to consider the coupled vibration of gear translation, bending, torsion, pendulum and other vibration forms. In order 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 Fig.11 As shown in the figure, local coordinate systems O1-X1Y1Z1 and O2-X2Y2Z2 are established on the driving and driven bevel gears respectively. Fig.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, and the meshing master and slave spiral bevel gears are simplified into two lumped mass points. The meshing relationship is represented by a spring connection between the gear pairs. 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 the gear transmission system, rolling bearings are usually used as the supporting 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 coupling dimensions 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] In the formula, 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.

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

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

[0132] Firstly, based on the gear tooth contact analysis model, the initial contact point and the coordinates of the tooth surface points of the large and small gears 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 the tooth surface deformation is extracted by applying a normal load on the contact surface tooth surface, and the flexibility matrix is ​​established. Then, the deformation coordination equation is established to solve the actual contact point and the tooth surface node distribution force.

[0133] Secondly, based on the tooth surface point coordinates of the large and small gears obtained by the gear contact analysis model, the geometric model of the spiral bevel gear is established through three-dimensional software; the geometric model is imported into the finite element software, and constraints are added and meshing is performed. After meshing, the tooth surface node coordinates are extracted, and the coordinates obtained from the gear contact analysis model are matched with the finite element tooth surface coordinates by finding the closest point. Based on the definition of ΔK in the Paris law, the variable tooth surface node distribution force obtained by the load-bearing gear contact analysis is simplified to the distributed force of fixed position cyclic loading. Based on the simplified distributed force and the tooth surface point coordinates of the finite element, the initial crack is inserted and the crack is extended, thereby establishing an efficient crack extension model. A single crack growth threshold is set. When the crack extends to the single crack growth threshold, the stiffness is calculated based on the gear model after crack extension.

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

[0135] The present invention has many specific application paths, and the above is only a preferred embodiment of the present invention. It should be pointed out that for ordinary technicians in this technical field, several improvements and modifications can be made without departing from the principle of the present invention, and these improvements and modifications should also be regarded as the protection scope 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 the gear tooth contact analysis; Potential contact points are calculated using contact ellipses; By applying normal load on the contact surface, the tooth surface deformation is extracted and the potential contact point flexibility matrix is ​​established; then, the deformation coordination equation is established to solve the tooth surface node distribution force; (2) The geometric model of the arc bevel gear pair is established through the tooth surface points of the large and small gears; the cyclically changing distributed force on one tooth surface of the arc bevel gear is simplified to the tooth surface point distributed force at a fixed meshing position; the crack propagation model is established through the geometric model of the arc bevel gear pair, the simplified tooth surface point distributed force at a fixed meshing position, and the initial crack is 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 extension in the next stage and obtain the shape of the crack after extension; A single crack growth threshold is set to perform crack extension. When the crack extends to the single crack growth threshold, the meshing stiffness of the finite element model of the spiral bevel gear pair containing cracks is calculated by the established load-bearing gear tooth contact analysis model; the finite element model of the spiral bevel gear pair containing cracks is a finite element model obtained after crack extension simulation based on the crack extension model; (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 introduced into the dynamic model; the meshing force of the tooth surface is calculated, and then the meshing force is introduced 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 re-substituted 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 through the crack propagation model, thereby forming a coupling process until the crack propagation ends.

2. The method for calculating the time-varying meshing stiffness of spiral bevel gears according to claim 1, characterized in that: In step (1), each initial contact point represents a meshing moment and a gear meshing position.

3. The method for calculating the time-varying meshing stiffness of spiral bevel gears according to claim 1, characterized in that: In step (1), a finite element model is established using a three-dimensional bilinear isoparametric unit; 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, and the deformation at a position below a certain distance from 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: In the formula, and are the corresponding element values ​​of the compliance 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; And calculate the contact compliance λ c : 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 calculating the time-varying meshing stiffness of spiral bevel gears according to claim 1, characterized in that: 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: In the formula, 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 meshing of the spiral bevel gear.

5. The method for calculating the time-varying meshing stiffness of spiral bevel gears according to claim 1, characterized in that: In step (2), a geometric model of a spiral bevel gear pair is established by three-dimensional software, and the geometric model of the spiral bevel gear pair is imported into a finite element software to set boundary conditions, divide the mesh, and extract the node coordinates of the finite element tooth surface; Simplify the force on the tooth surface based on the Paris formula; 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, program through MATLAB software 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, match the potential contact point with the finite element tooth surface node, and apply the simplified normal distribution force to the finite element tooth surface node through parametric programming; Introduce the initial crack into the root of the gear tooth and perform meshing simulation, calculate the stress intensity factor distribution at the crack tip based on the stress intensity factor theory, calculate the extension step based on the Paris formula, and calculate the extension direction based on the MTS theory, so as to perform crack extension; Substitute the finite element model after crack propagation into the finite element software and recalculate the meshing to obtain the transmission error.

6. The method for calculating the time-varying meshing stiffness of spiral bevel gears according to claim 5, characterized in that: 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.

7. The method for calculating the time-varying meshing stiffness of spiral bevel gears according to claim 6, characterized in that: The meshing stiffness K calculated based on the above crack extension is substituted into the dynamic model to solve and obtain the required tooth surface meshing force.

8. The method for calculating the time-varying meshing stiffness of spiral bevel gears according to claim 7, characterized in that: According to the stiffness matrix assembly method in the finite element method, various units are assembled into the overall system matrix in sequence according to the unit node numbers shown in the figure above. After integrating all matrices, the dynamic model of the entire system is obtained: In the formula, 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 second derivative of X is the acceleration vector.

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

10. 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 8 are implemented.

Images