Tooth surface friction calculation method based on mathematical programming algorithm

Abstract

The invention belongs to the field of machinery, and particularly relates to a tooth surface friction calculation method based on a mathematical programming algorithm. The objective of the invention is to solve the problem of low calculation efficiency of the existing numerical analysis technology for the gear contact problem. According to the tooth surface friction calculation method based on the mathematical programming algorithm, the problem scale is greatly reduced through substructure static condensation, and the high efficiency of the mathematical programming algorithm (the smooth Newton method) is combined, so that the leap of the calculation speed is achieved. In addition, according to the method, a mechanical problem is creatively converted into a mathematical programming problem, and a brand new technical path is provided. The method is not only suitable for straight teeth, helical teeth and herringbone gears, but also can be popularized to friction contact analysis of other mechanical connection pairs (such as splines and bearing raceways) in theoretical framework. The whole process from parametric modeling, data exporting to model solving can achieve high automation through scripts, human intervention is reduced, and the repeatability of results is guaranteed.

Description

Technical Field

[0001] This invention belongs to the field of mechanics, and in particular relates to a method for calculating tooth surface friction based on mathematical programming algorithms. Background Technology

[0002] Gear transmission is one of the most important forms of power and motion transmission in modern mechanical equipment. The contact stress and friction behavior between the gear teeth directly affect the gear's load-bearing capacity, transmission efficiency, vibration noise, and service life. Therefore, accurate numerical simulation of the tooth surface friction contact problem is crucial for gear design and performance evaluation.

[0003] Currently, general-purpose finite element software (such as ANSYS and Abaqus) is widely used in engineering to perform numerical analysis of gear contact problems. These software programs typically provide the following two main types of contact algorithms:

[0004] 1. Penalty Function Method: This method simulates contact by introducing a virtual spring, where the contact force is proportional to the penetration amount. Its solution accuracy heavily depends on the user-defined "penalty factor." A penalty factor that is too small leads to excessive physical penetration and distorted results; a penalty factor that is too large deteriorates the behavior of the system equations, causing convergence difficulties. Therefore, this method can only approximately satisfy the contact conditions, and its accuracy and stability are difficult to guarantee.

[0005] 2. Lagrange multiplier method: This method introduces additional Lagrange multipliers (i.e., contact forces) as degrees of freedom for the system, which can precisely satisfy the impenetrability condition on the contact surface. However, this leads to an increase in the total number of degrees of freedom of the system. Especially when dealing with complex three-dimensional contact problems involving friction, the computational scale expands dramatically, the solution time is very long, and the computational cost is high.

[0006] The shortcomings of existing numerical analysis techniques for gear contact problems:

[0007] (1) Low computational efficiency: Traditional finite element contact analysis, especially three-dimensional problems with friction, requires a large number of iterative calculations, which is time-consuming and difficult to meet the needs of modern rapid design and optimization iteration.

[0008] (2) It is difficult to balance accuracy and stability: the accuracy of the penalty function method is greatly affected by human parameters; although the Lagrange multiplier method has high accuracy, it has a high computational cost.

[0009] (3) Poor adaptability to complex contact boundaries: During gear meshing, the contact area changes constantly, and the search and judgment of contact pairs requires a lot of computing resources.

[0010] Therefore, existing numerical analysis techniques for gear contact problems suffer from low computational efficiency. Summary of the Invention

[0011] The purpose of this invention is to address the low computational efficiency of existing numerical analysis techniques for gear contact problems. It provides a method for calculating tooth surface friction based on a mathematical programming algorithm, comprising:

[0012] S1: Obtain the geometric parameter data of the driving gear and the driven gear to be calculated;

[0013] Based on the geometric parameter data of the driving gear and the driven gear to be calculated, mathematical calculation software is used to calculate the tooth profile node coordinate data of the driving gear and the tooth profile node coordinate data of the driven gear.

[0014] S2: Based on the tooth profile node coordinate data of the driving gear and the tooth profile node coordinate data of the driven gear, parametric modeling is performed using finite element software to obtain the finite element model of the driving gear and the finite element model of the driving gear.

[0015] S3: In the finite element software, the finite element models of the driving gear and the driving gear are assembled according to the meshing principle to obtain the gear pair meshing finite element model.

[0016] S4: Determine the potential contact area of ​​the gear pair meshing finite element model, refine the potential contact area using the grid method in the finite element software, and define the nodes of the refined contact area as the master degree of freedom nodes of the tooth surface of the gear pair meshing finite element model.

[0017] S5: In finite element software, the substructure static condensation method is used to condense all degrees of freedom in the gear pair meshing finite element model onto the main degree of freedom nodes on the tooth surface to obtain a super element model.

[0018] The super-element model includes: a driving wheel super-element model and a driven wheel super-element model;

[0019] The substructure static condensation method is achieved through the super-element function in finite element software. It condenses the degrees of freedom of the driving gear in the gear pair meshing finite element model to the main degree-of-freedom nodes of the driving gear tooth surface, and condenses the degrees of freedom of the driven gear super-element in the gear pair meshing finite element model to the main degree-of-freedom nodes of the driven gear super-element tooth surface. This is a well-known process among those skilled in the art.

[0020] Export the stiffness matrix of the driving wheel, the first load vector of the driving wheel, the second load vector of the driving wheel, and the nodal coordinate information of the driving wheel in the driving wheel super-element model;

[0021] Export the driven wheel stiffness matrix, the first load vector of the driven wheel, the second load vector of the driven wheel, and the nodal coordinate information of the driven wheel in the driven wheel super-element model;

[0022] The specific process is as follows:

[0023] Using non-public commands from finite element software (such as ANSYS's / DEBUG command), the condensed stiffness matrix, load vectors, and nodal coordinate information were exported as text files. Subsequently, a Python script was used for data cleaning and format conversion to generate a .mat file readable by MATLAB.

[0024] S6: Based on the stiffness matrix of the driving wheel, the first load vector of the driving wheel, the second load vector of the driving wheel, the nodal coordinate information of the driving wheel, the stiffness matrix of the driven wheel, the first load vector of the driven wheel, the second load vector of the driven wheel, and the nodal coordinate information of the driven wheel, construct a point-plane second-order conical linear complementary model;

[0025] The second-order conical linear complementary model uses generalized contact force and generalized relative displacement as complementary variables;

[0026] S7: The second-order conical linear complementary model is solved using the smooth Newton algorithm in mathematical calculation software to obtain the tooth surface friction calculation results, including: the normal contact force, tangential friction force, normal clearance and tangential slip of the tooth surface;

[0027] The gear geometric parameter data to be calculated in S1 includes: involute data, addendum circle data, and dedendum circle data.

[0028] Preferably, the mathematical calculation software used in this invention can be MATLAB;

[0029] In step S1, based on the geometric parameter data of the driving gear and the driven gear to be calculated, mathematical calculation software is used to calculate the tooth profile node coordinate data of the driving gear and the driven gear. The specific process is as follows:

[0030] Input the geometric parameters of the gear to be calculated and the geometric parameters of the driven gear to be calculated into the mathematical calculation software;

[0031] Based on the geometric parameters of the driving gear to be calculated, the parameterized modeling function of mathematical calculation software is used to discretely generate the tooth profile node coordinate data of the driving gear.

[0032] Based on the geometric parameter data of the driven gear to be calculated, the parameterized modeling function of mathematical calculation software is used to discretely generate the tooth profile node coordinate data of the driven gear.

[0033] In step S2, based on the tooth profile node coordinate data of the driving gear and the driven gear, parametric modeling is performed using finite element software to obtain the finite element model of the driving gear and the finite element model of the driving gear; the specific process is as follows:

[0034] Input the tooth profile node coordinate data of the driving gear and the driven gear obtained from S1 into the finite element software;

[0035] Based on the tooth profile node coordinate data of the driving gear, the node-element direct generation function in the finite element software is used to generate the finite element model of the driving gear. The specific process is as follows: the tooth profile node coordinate data of the driving gear are sequentially connected using the node connection method to generate tooth profile end face elements. Then, the VDRAG extrusion command is used to complete the establishment of the solid element in the tooth width direction of a single tooth. Finally, the VGEN rotation complex command is used on the single tooth element to complete the establishment of the finite element model of the driving gear.

[0036] Based on the tooth profile node coordinate data of the driven gear, the node-element direct generation function in the finite element software is used to generate the finite element model of the driven gear. The specific process is as follows: the node connection method is used to sequentially connect the tooth profile node coordinate data of the driven gear to generate tooth profile end face elements. Then, the VDRAG extrusion command is used to complete the establishment of the solid element in the tooth width direction of a single tooth. Finally, the VGEN rotation complex command is used on the single tooth element to complete the establishment of the finite element model of the driven gear.

[0037] In step S3, the finite element models of the driving gear and the driven gear are assembled according to the meshing principle to obtain the gear pair meshing finite element model. The specific process is as follows:

[0038] S3.1: Based on the gear meshing principle, set the initial meshing phase of the driving gear and the initial meshing phase of the driven gear;

[0039] S3.2: Set the constraint conditions for the finite element model of the driven gear. The specific process is as follows: Apply full constraints to the driven gear.

[0040] S3.3: Set the constraints for the finite element model of the driving gear. The specific process is as follows: constrain the driving gear for all degrees of freedom except for the axial rotational degree of freedom, and apply the working torque; to obtain the gear pair meshing finite element model.

[0041] In step S6, a point-plane second-order conical linear complementary model is constructed based on the stiffness matrix of the driving wheel, the first load vector of the driving wheel, the second load vector of the driving wheel, the node coordinate information of the driving wheel, the stiffness matrix of the driven wheel, the first load vector of the driven wheel, the second load vector of the driven wheel, and the node coordinate information of the driven wheel; the specific process is as follows:

[0042] S6.1: Based on the coordinate information of the driving wheel node and the driven wheel node, K point-to-surface friction contact pairs are constructed in the mathematical calculation software. Each point-to-surface friction contact pair consists of a contact node and a target surface.

[0043] The contact node is a node on the active gear tooth surface, and the target surface is a four-node isoparametric element on the driven gear tooth surface;

[0044] The four-node isoparametric element is a standardized mathematical tool used to precisely define and describe the geometry and deformation behavior of a small region (i.e., the "target surface") on the tooth surface of the driven gear. It can be understood as a virtual patch with specific rules that can be flexibly deformed. This definition allows for contact mesh mismatch; a four-node isoparametric element consists of four nodes forming a surface patch; a four-node isoparametric element is composed of four points (K1, K2, K3, K4) on the tooth surface of the driven gear. In the actual gear finite element model, these four points are four adjacent nodes on the tooth surface mesh. These four nodes uniquely determine a surface patch (usually a slightly curved quadrilateral) in three-dimensional space, which represents a small potential contact area on the gear tooth surface.

[0045] Point-to-surface contact elements are used to describe the tooth surface friction contact model. This allows for the simulation of frictional contact under unloaded conditions. Such planar elements can effectively represent the potential contact area of ​​the target surface. On the other hand, it is assumed that different nodes on the contact surface (such as...) (Point) contacts the target surface. Dot at On a plane, and Point to plane The shortest distance. In this way, the contact area mesh does not need to be matched into a one-to-one correspondence between nodes, and the contact constraint relationship between the two tooth surfaces can be established through point-to-surface contact elements on the contact system.

[0046] Due to the position of the point of action on the target surface It is unknown. As we know from the theoretical foundation in Chapter 2, certain assumptions need to be set for the model:

[0047] 1) The contact area between the two tooth surfaces is within the linear elastic range, which conforms to the small deformation theory in elasticity mechanics, and kinetic friction is negligible;

[0048] 2) The gap between the two tooth surfaces at the possible contact boundary is small enough. Point and The spatial distance between points is And satisfy yes Perpendicular to the target surface:

[0049] 3) The frictional contact between the two tooth surfaces satisfies the classical Coulomb friction law.

[0050] All contact pairs in the contact system are assumed to satisfy the above conditions, and the contact pairs satisfy the material's normal impermeability condition, while the contact friction and tangential slip conform to Coulomb's law of friction. To ensure that the contact system consists of reasonable point-to-surface contact pairs, multiple contact points can simultaneously form different contact pairs with a target surface.

[0051] S6.2: Establish a local coordinate system for each point-to-surface friction contact pair element to obtain K local coordinate systems;

[0052] And define a generalized contact force vector and a generalized relative displacement vector in each local coordinate system;

[0053] S6.3: Construct the transformation matrix based on the K local coordinate systems;

[0054] S6.4: Based on the transformation matrix, the stiffness matrix of the driving wheel, the first load vector of the driving wheel, and the second load vector of the driving wheel, construct the equilibrium equation of the driving wheel in the local coordinate system;

[0055] S6.5: Based on the transformation matrix, the driven wheel stiffness matrix, the driven wheel first load vector, and the driven wheel second load vector, construct the driven wheel equilibrium equation in the local coordinate system;

[0056] S6.6: Construct a second-order conical linear complementary model by combining the equilibrium equations of the driving wheel in the local coordinate system and the equilibrium equations of the driven wheel in the local coordinate system.

[0057] The beneficial effects of this invention are as follows:

[0058] This invention significantly reduces the problem size through substructure static condensation, and combined with the efficiency of mathematical programming algorithms (smooth Newton's method), achieves a leap in computational speed. Experimental results show that the total computation time is only 1 / 4 of that of traditional ANSYS contact analysis (an improvement of approximately 4.22 times), making it particularly suitable for parametric studies and optimization design. Based on rigorous mathematical theory, the contact friction condition is equivalent to a second-order cone complementary condition, thus accurately satisfying Coulomb's friction law and the normal contact condition, avoiding the approximation errors of the penalty function method. Compared with commercial software, the maximum relative error is less than 10%, and the agreement is even higher in non-stress concentration regions.

[0059] Furthermore, this invention creatively transforms a mechanical problem into a mathematical programming problem, providing a completely new technical approach. This method is not only applicable to spur, helical, and herringbone gears, but its theoretical framework can also be extended to the frictional contact analysis of other mechanical connection pairs (such as splines and bearing raceways). From parametric modeling and data export to model solving, the entire process can be highly automated through scripts, reducing human intervention and ensuring the repeatability of results. Attached Figure Description

[0060] Figure 1 This is a schematic model of the tooth surface friction contact problem of the present invention;

[0061] Figure 2 This is a schematic diagram of the contact unit and local coordinate system of the present invention;

[0062] Figure 3 This is a schematic diagram of the 4-node isoparametric element of the target surface of the present invention;

[0063] Figure 4 This is a schematic diagram of the solution process for the second-order conical linear complementary model of this invention;

[0064] Figure 5 This is a schematic diagram of the maximum frictional stress distribution in the tooth width direction obtained by solving the second-order cone model of this invention;

[0065] Figure 6 This is a schematic diagram of the projection of the contact stress on the left-hand tooth surface of the present invention onto the meshing plane;

[0066] Figure 7 This is a schematic diagram of the projection of the contact stress on the right-hand tooth surface of the present invention onto the meshing plane;

[0067] Figure 8 This is a schematic diagram of the projection of the frictional stress on the left-hand tooth surface of the present invention onto the meshing plane;

[0068] Figure 9 This is a schematic diagram of the projection of the frictional stress on the right-hand tooth surface of the present invention onto the meshing plane. Detailed Implementation

[0069] Specific implementation method one: Combining Figures 1-9 This invention is described;

[0070] S1: Obtain the geometric parameter data of the driving gear and the driven gear to be calculated;

[0071] Based on the geometric parameter data of the driving gear and the driven gear to be calculated, mathematical calculation software is used to calculate the tooth profile node coordinate data of the driving gear and the tooth profile node coordinate data of the driven gear.

[0072] S2: Based on the tooth profile node coordinate data of the driving gear and the tooth profile node coordinate data of the driven gear, parametric modeling is performed using finite element software to obtain the finite element model of the driving gear and the finite element model of the driving gear.

[0073] S3: In the finite element software, the finite element models of the driving gear and the driving gear are assembled according to the meshing principle to obtain the gear pair meshing finite element model.

[0074] S4: Determine the potential contact area of ​​the gear pair meshing finite element model, refine the potential contact area using the grid method in the finite element software, and define the nodes of the refined contact area as the master degree of freedom nodes of the tooth surface of the gear pair meshing finite element model.

[0075] S5: In finite element software, the substructure static condensation method is used to condense all degrees of freedom in the gear pair meshing finite element model onto the main degree of freedom nodes on the tooth surface to obtain a super element model.

[0076] The super-element model includes: a driving wheel super-element model and a driven wheel super-element model;

[0077] The substructure static condensation method is achieved through the super-element function in finite element software. It condenses the degrees of freedom of the driving gear in the gear pair meshing finite element model to the main degree-of-freedom nodes of the driving gear tooth surface, and condenses the degrees of freedom of the driven gear super-element in the gear pair meshing finite element model to the main degree-of-freedom nodes of the driven gear super-element tooth surface. This is a well-known process among those skilled in the art.

[0078] Export the stiffness matrix of the driving wheel, the first load vector of the driving wheel, the second load vector of the driving wheel, and the nodal coordinate information of the driving wheel in the driving wheel super-element model;

[0079] Export the driven wheel stiffness matrix, the first load vector of the driven wheel, the second load vector of the driven wheel, and the nodal coordinate information of the driven wheel in the driven wheel super-element model;

[0080] The specific process is as follows:

[0081] Using non-public commands from finite element software (such as ANSYS's / DEBUG command), the condensed stiffness matrix, load vectors, and nodal coordinate information were exported as text files. Subsequently, a Python script was used for data cleaning and format conversion to generate a .mat file readable by MATLAB.

[0082] S6: Based on the stiffness matrix of the driving wheel, the first load vector of the driving wheel, the second load vector of the driving wheel, the nodal coordinate information of the driving wheel, the stiffness matrix of the driven wheel, the first load vector of the driven wheel, the second load vector of the driven wheel, and the nodal coordinate information of the driven wheel, construct a point-plane second-order conical linear complementary model;

[0083] The second-order conical linear complementary model uses generalized contact force and generalized relative displacement as complementary variables;

[0084] S7: The second-order conical linear complementary model is solved using the smooth Newton algorithm in mathematical calculation software to obtain the tooth surface friction calculation results, including: the normal contact force, tangential friction force, normal clearance, and tangential slip of the tooth surface.

[0085] Specific Implementation Method Two: The difference between this implementation method and Specific Implementation Method One is that:

[0086] The gear geometric parameter data to be calculated in S1 includes: involute data, addendum circle data, and dedendum circle data.

[0087] Preferably, the mathematical calculation software used in this invention can be MATLAB;

[0088] In step S1, based on the geometric parameter data of the driving gear and the driven gear to be calculated, mathematical calculation software is used to calculate the tooth profile node coordinate data of the driving gear and the driven gear. The specific process is as follows:

[0089] Input the geometric parameters of the gear to be calculated and the geometric parameters of the driven gear to be calculated into the mathematical calculation software;

[0090] Based on the geometric parameters of the driving gear to be calculated, the parameterized modeling function of mathematical calculation software is used to discretely generate the tooth profile node coordinate data of the driving gear.

[0091] Based on the geometric parameters of the driven gear to be calculated, the parametric modeling function of mathematical calculation software is used to discretize and generate the tooth profile node coordinate data of the driven gear.

[0092] The other steps and parameters are the same as in Specific Implementation Method 1.

[0093] Specific Implementation Method Three: The difference between this implementation method and Specific Implementation Methods One and Two is that:

[0094] In step S2, based on the tooth profile node coordinate data of the driving gear and the driven gear, parametric modeling is performed using finite element software to obtain the finite element model of the driving gear and the finite element model of the driving gear; the specific process is as follows:

[0095] Input the tooth profile node coordinate data of the driving gear and the driven gear obtained from S1 into the finite element software;

[0096] Based on the tooth profile node coordinate data of the driving gear, the node-element direct generation function in the finite element software is used to generate the finite element model of the driving gear. The specific process is as follows: the tooth profile node coordinate data of the driving gear are sequentially connected using the node connection method to generate tooth profile end face elements. Then, the VDRAG extrusion command is used to complete the establishment of the solid element in the tooth width direction of a single tooth. Finally, the VGEN rotation complex command is used on the single tooth element to complete the establishment of the finite element model of the driving gear.

[0097] Based on the tooth profile node coordinate data of the driven gear, the node-element direct generation function in the finite element software is used to generate the finite element model of the driven gear. The specific process is as follows: the node connection method is used to sequentially connect the tooth profile node coordinate data of the driven gear to generate tooth profile end face elements. Then, the VDRAG extrusion command is used to complete the establishment of the solid element in the tooth width direction of a single tooth. Finally, the VGEN rotation complex command is used on the single tooth element to complete the establishment of the finite element model of the driven gear.

[0098] The other steps and parameters are the same as in one of the specific implementation methods one or two.

[0099] Specific Implementation Method Four: This implementation method differs from Specific Implementation Methods One to Three in that:

[0100] In step S3, the finite element models of the driving gear and the driven gear are assembled according to the meshing principle to obtain the gear pair meshing finite element model. The specific process is as follows:

[0101] S3.1: Based on the gear meshing principle, set the initial meshing phase of the driving gear and the initial meshing phase of the driven gear;

[0102] S3.2: Set the constraint conditions for the finite element model of the driven gear. The specific process is as follows: Apply full constraints to the driven gear.

[0103] S3.3: Set the constraint conditions for the finite element model of the driving gear. The specific process is as follows: constrain the driving gear except for the axial rotational degree of freedom, and apply the working torque; to obtain the gear pair meshing finite element model.

[0104] The other steps and parameters are the same as those in one of the specific implementation methods one to three.

[0105] Specific Implementation Method Five: The difference between this implementation method and Specific Implementation Methods One to Four is that:

[0106] In step S6, a point-plane second-order conical linear complementary model is constructed based on the stiffness matrix of the driving wheel, the first load vector of the driving wheel, the second load vector of the driving wheel, the node coordinate information of the driving wheel, the stiffness matrix of the driven wheel, the first load vector of the driven wheel, the second load vector of the driven wheel, and the node coordinate information of the driven wheel; the specific process is as follows:

[0107] S6.1: Based on the coordinate information of the driving wheel node and the driven wheel node, K point-to-surface friction contact pairs are constructed in the mathematical calculation software. Each point-to-surface friction contact pair consists of a contact node and a target surface.

[0108] The contact node is a node on the active gear tooth surface, and the target surface is a four-node isoparametric element on the driven gear tooth surface;

[0109] The four-node isoparametric element is a standardized mathematical tool used to precisely define and describe the geometry and deformation behavior of a small region (i.e., the "target surface") on the tooth surface of the driven gear. It can be understood as a virtual patch with specific rules that can be flexibly deformed. This definition allows for contact mesh mismatch; a four-node isoparametric element consists of four nodes forming a surface patch; a four-node isoparametric element is composed of four points (K1, K2, K3, K4) on the tooth surface of the driven gear. In the actual gear finite element model, these four points are four adjacent nodes on the tooth surface mesh. These four nodes uniquely determine a surface patch (usually a slightly curved quadrilateral) in three-dimensional space, which represents a small potential contact area on the gear tooth surface.

[0110] Point-to-surface contact elements are used to describe the tooth surface friction contact model. This allows for the simulation of frictional contact under unloaded conditions. Such planar elements can effectively represent the potential contact area of ​​the target surface. On the other hand, it is assumed that different nodes on the contact surface... (like (Point) contacts the target surface. Dot at On a plane, and Point to plane The shortest distance. In this way, the contact area mesh does not need to be matched into a one-to-one correspondence between nodes, and the contact constraint relationship between the two tooth surfaces can be established through point-to-surface contact elements on the contact system.

[0111] Due to the position of the point of action on the target surface It is unknown. As we know from the theoretical foundation in Chapter 2, certain assumptions need to be set for the model:

[0112] (1) The contact body of the two tooth surfaces is in the linear elastic range, which conforms to the small deformation theory in elasticity, and kinetic friction is ignored;

[0113] (2) The gap between the two tooth surfaces at the possible contact boundary is small enough. Point and The spatial distance between points is And satisfy yes Perpendicular to the target surface:

[0114] (3) The frictional contact between the two tooth surfaces satisfies the classical Coulomb friction law.

[0115] All contact pairs in the contact system are assumed to satisfy the above conditions, and the contact pairs satisfy the material's normal impermeability condition, while the contact friction and tangential slip conform to Coulomb's law of friction. To ensure that the contact system consists of reasonable point-to-surface contact pairs, multiple contact points can simultaneously form different contact pairs with a target surface.

[0116] S6.2: Establish a local coordinate system for each point-to-surface friction contact pair element to obtain K local coordinate systems;

[0117] And define a generalized contact force vector and a generalized relative displacement vector in each local coordinate system;

[0118] S6.3: Construct the transformation matrix based on the K local coordinate systems;

[0119] S6.4: Based on the transformation matrix, the stiffness matrix of the driving wheel, the first load vector of the driving wheel, and the second load vector of the driving wheel, construct the equilibrium equation of the driving wheel in the local coordinate system;

[0120] S6.5: Based on the transformation matrix, the driven wheel stiffness matrix, the driven wheel first load vector, and the driven wheel second load vector, construct the driven wheel equilibrium equation in the local coordinate system;

[0121] S6.6: Construct a second-order conical linear complementary model by combining the equilibrium equations of the driving wheel and the driven wheel in the local coordinate system;

[0122] The other steps and parameters are the same as those in one of the specific implementation methods one to four.

[0123] Specific Implementation Method Six: The difference between this implementation method and Specific Implementation Methods One through Five is that:

[0124] In step S6.2, a local coordinate system is established for each point-to-surface friction contact pair, resulting in K local coordinate systems; the specific process is as follows:

[0125] To establish a local coordinate system for the k-th point-to-surface friction contact pair element, the specific process for obtaining the k-th local coordinate system is as follows:

[0126] After the driving wheel is subjected to working torque, the contact surface node and the target surface of the driven wheel come into frictional contact and produce elastic deformation.

[0127] Let the contact nodes in the k-th point-to-surface friction contact element be at a certain direction angle. Along Move in the direction normal, and in the direction relative to the target surface. and Tangential slip, , and Both are unit vectors.

[0128] According to unit vector Unit vector With unit vector Construct the k-th local spatial coordinate system;

[0129] like Figure 3 As shown, the elastic deformation of each contact node is converted into a local normal vector and a tangent vector, which facilitates the establishment of the subsequent frictional contact state equation; each frictional contact pair element satisfies the following contact control equation:

[0130] (1)

[0131] It is the generalized slip vector of the k-th friction contact pair element, where It is the total tangential slip of the contact pair element. It is the tangential slip of the contact pair element in the two tangential directions of the local coordinate system. It is the normal gap between contact pairs; It is the generalized contact force vector of the k-th frictional contact pair element, where It is the normal contact force between the contact elements. It is the tangential frictional force between the contact pair elements in two tangential directions in the local coordinate system. It is the coefficient of friction between the two tooth surfaces; Let represent the Cartesian product of a three-dimensional second-order cone and a one-dimensional non-negative cone.

[0132] According to the finite element theory, the equilibrium equations of the two-gear contact system can be discretized and expressed as follows: .

[0133] in: The total stiffness matrix; The nodal displacement vectors are discrete representations of the model. This is the node equivalent load vector.

[0134] The finite element equations above are expanded, rearranged, and sorted, dividing all degrees of freedom into potential contact node degrees of freedom and other non-possible contact node degrees of freedom. The corresponding total stiffness matrix and nodal equivalent load vectors are also rearranged and divided into blocks:

[0135] (2)

[0136] Where: subscript 1 represents the potential contact surface node of the driving wheel, subscript 2 represents the potential target surface node of the driven wheel, and subscript 3 represents the remaining non-possible contact nodes of the model.

[0137] The first static condensation is performed using the substructure method in the finite element method, which condenses the finite element model of the model onto the potential contact boundary of the tooth surface. The potential contact boundary of the entire model only occupies a small part of the system. Static condensation can reduce the computational dimension of the entire model.

[0138] First, the displacement vectors of the remaining non-contact nodes are expressed as: .

[0139] Then, substituting the above equation into the balance equations of the first two dimensions of the model and rearranging, we obtain:

[0140]

[0141] in: , , ,

[0142] , , ,

[0143] Since the potential contact surface node of the driving gear and the potential target surface node of the driven gear are two independent elastic gear entities, the stiffness matrix regarding their interaction relationship is... and It is actually a zero matrix;

[0144] At the same time, the equivalent load vector of the nodes is reorganized and decomposed. Decomposed into the first load vector (a generalized vector representing the external loads of the nodes). ) and the second load vector (a generalized vector representing the contact force between the contact pairs at the nodes) Therefore, the equilibrium equation after static condensation reduction is:

[0145] (4)

[0146] The above-mentioned substructure method for static condensation operations can be implemented in ANSYS by defining superelements. After performing substructure calculations, the / DEBUG command can be used.

[0147] The substructure stiffness matrix of the potential contact node on the tooth surface includes: (Represents the stiffness matrix of the drive wheel) (Represents the stiffness matrix of the driven wheel)

[0148] The load array includes: (Represents the first load vector of the drive wheel) (Represents the first load vector of the driven wheel) (Second load vector of the drive wheel) (Second load vector of driven wheel);

[0149] Export the above data as a TXT text file, then use a Python program to clean and convert the data to a Mat data format that can be directly read by Matlab.

[0150] The other steps and parameters are the same as those in any of the specific implementation methods one to five.

[0151] Specific Implementation Method Seven: The difference between this implementation method and Specific Implementation Methods One through Six is ​​that:

[0152] In S6.3, the transformation matrix is ​​constructed based on K local coordinate systems, which can be expressed by the following formula:

[0153]

[0154] In the formula, Represents the transformation matrix. express The element in the first row and first column, express The element in the k-th row and k-th column, express The element in the Kth row and Kth column, the The element in row k and column k The calculation formula is expressed as:

[0155]

[0156] In the formula, x represents the x-direction vector, y represents the y-direction vector, and z represents the z-direction vector. Represents the cosine value

[0157] The other steps and parameters are the same as those in any of the specific implementation methods one to six.

[0158] Specific Implementation Method Eight: The difference between this implementation method and Specific Implementation Methods One through Seven is that:

[0159] In step S6.4, the equilibrium equations of the driving wheel in the local coordinate system are constructed based on the transformation matrix, the stiffness matrix of the driving wheel, the first load vector of the driving wheel, and the second load vector of the driving wheel; expressed by the formula:

[0160]

[0161] In the formula, This represents the stiffness matrix of the driving wheel in the local coordinate system. This represents the first load vector of the driving wheel in the local coordinate system. This represents the second load vector of the driving wheel in the local coordinate system. The displacement vector of the driving wheel node in the local coordinate system is expressed by the formula:

[0162]

[0163]

[0164]

[0165]

[0166] In the formula, express The inverse matrix, Represents the stiffness matrix of the drive wheel, Indicates the first load vector of the driving wheel, This represents the second load vector of the driving wheel; Represents the displacement vector of the active wheel node;

[0167] (a) The derivation process of the equilibrium equation of the driving wheel in the local coordinate system is as follows: Expanding the equilibrium equation after static condensation, the equilibrium equation of the potential contact node of the driving wheel in the global coordinate system is:

[0168] (5)

[0169] In the formula: It is the stiffness matrix of the contact surface; as well as It is a generalized column vector of external loads and contact forces acting on the contact surface; Let the displacement vector of the nodes on the potential contact surface of the driving wheel be the following equation:

[0170] (6)

[0171] In the formula: Represented as the contact surface, Both K and K represent the number of potential contact nodes on the driving wheel. K is used in a general sense. The two variables are used to distinguish between the number of nodes of the driven wheel and the number of nodes of the driven wheel; the two variables have the same meaning.

[0172] After the driving wheel is subjected to the working torque, the contact surface nodes of the driving wheel and the target surface of the driven wheel undergo elastic deformation due to frictional contact. Assume the k-th contact node is at a certain direction angle. Along Move in the direction normal, and in the direction relative to the target surface. and Tangential slip, and It is a unit vector, and is related to the unit vector. A local spatial coordinate system is constructed. The elastic deformation of each contact node is converted into a local normal vector and a tangent vector, which facilitates the subsequent establishment of the friction contact state equation. Therefore, the transformation matrix between the local coordinate system and the global coordinate system is utilized. Transform the nodal degrees of freedom of the contact surface to the local contact coordinate system:

[0173] (7)

[0174] In the formula: ;

[0175]

[0176] .

[0177] Similarly, the generalized column vectors of external loads and contact forces acting on the contact surface are transformed into the local contact coordinate system:

[0178] (8)

[0179] in: ;

[0180] ;

[0181] ;

[0182] .

[0183] Substituting the generalized displacement vector, external load vector, and contact force vector, transformed into the contact pair element coordinate system, into the equilibrium equation of the potential contact surface of the driving wheel, we obtain the following equation:

[0184] (9)

[0185] The transformation matrix is ​​an orthogonal matrix; multiplying both sides of the above equation by... The equilibrium equations of the potential contact surface of the driving wheel are then expressed in local coordinates:

[0186] (10)

[0187] Wherein: the stiffness matrix of the potential contact surface of the driving wheel in the local coordinate system is:

[0188] .

[0189] The other steps and parameters are the same as those in any of the specific implementation methods one to seven.

[0190] Specific Implementation Method Nine: The difference between this implementation method and Specific Implementation Methods One through Eight is that:

[0191] In S6.5, based on the transformation matrix, the driven wheel stiffness matrix, the driven wheel first load vector, and the driven wheel second load vector, the driven wheel equilibrium equation in the local coordinate system is constructed, which is expressed by the following formula:

[0192]

[0193] In the formula, This represents the stiffness matrix of the driven wheel in the local coordinate system. This represents the displacement vector of the driven wheel node in the local coordinate system. The second load vector of the driven wheel in the local coordinate system is expressed by the formula:

[0194]

[0195]

[0196]

[0197] In the formula, This represents the stiffness matrix of the driven wheel. This represents the displacement vector of the driven wheel node. This represents the second load vector of the driven wheel. Represents a pointer matrix. Represents a shape function matrix,

[0198] pointer matrix Expressed as a formula:

[0199]

[0200] In the formula, express The first element, express The kth element, express The Kth element,

[0201] Shape function matrix Expressed as a formula:

[0202]

[0203] express, express The element in the first row and first column, express The element in the k-th row and k-th column, express The element in row K and column K,

[0204] (b) The derivation process of constructing the equilibrium equation of the driven wheel in the local coordinate system is as follows: Expanding the equilibrium equation after static condensation, the equilibrium equation of the potential target surface of the driven wheel in the global coordinate system is:

[0205] (11)

[0206] in: It is the stiffness matrix of the target surface; as well as It is a generalized column vector of the external loads and contact forces acting on the target surface; For the nodal displacement vector on the potential target surface of the driven wheel, it can be extended to the following equation:

[0207] (12)

[0208] In the formula: Represented as the target surface, It is the number of nodes on the potential target surface of the driven wheel.

[0209] One contact node in the potential contact surface of the drive wheel A target surface on the potential target of the driven wheel, which is opposite to it. Defined as point-to-surface contact pair elements in the model.

[0210] Each 4-node target surface is indicated by a pointer matrix to show the position of the 4 nodes on all potential nodes of the driven wheel. The specific relationship is as follows:

[0211] (13)

[0212] in: The corresponding k is the number of point-to-surface contact pairs. .

[0213] The pointer matrix corresponding to each 4-node unit is as follows, and the dimension of the pointer matrix for each node is... ,

[0214] If this node is If a node is at the m-th position, then the corresponding pointer matrix is ​​at the m-th position. Listed to number The columns are identity matrices, and the remaining elements are 0, as shown below.

[0215] (14)

[0216] All target surface 4-node elements can be represented by a pointer matrix and all potential nodes of the driven wheel:

[0217] (15)

[0218] In the model, each 4-node target surface element belongs to a 4-node isoparametric element, such as... Figure 3 As shown, according to the isoparametric element theory in the finite element method, any quadrilateral element can be mapped to a square parametric element with a side length of 2, and a reference coordinate system can be established on it. In physical space, they do not need to be orthogonal to each other, nor do they need to be parallel to the global coordinate system; the edges of the unit are... and They bisect each other, and their equations are respectively and .

[0219] Target surface element With the target point The coordinates of the four nodes can be obtained in the global coordinate system. The global coordinates of the target point are Substituting the coordinates into the following formula yields the relative position of the target point within the square parameters. Then the shape function is calculated. The specific process is as follows:

[0220] (16)

[0221] Where: the shape function is The formula for calculating each element in the shape function is:

[0222] (17)

[0223] The formula for the coordinate vector of the corresponding 4-node isoparametric element is:

[0224] (18)

[0225] The frictional contact between the two tooth surfaces produces elastic deformation, target point respectively along The direction has shifted. This can be determined based on the nodal displacement. , and The displacement of the contact point is obtained. The following equation is satisfied:

[0226] (19)

[0227] in: .

[0228] Target surface element With the target point The relationship between the nodal displacement vectors is linked by shape functions, and can be expressed by the formula:

[0229] (20)

[0230] in: , It is a 3×3 identity matrix.

[0231] The nodal displacement vector corresponding to each target point is .

[0232] The target point is obtained by combining equations (13) and (20). The relationship between the global nodes of the potential target surface of the driven wheel and the global nodes is expressed by the formula:

[0233] (twenty one)

[0234] Similarly, to facilitate the establishment of the subsequent frictional contact state equations, the transformation matrix between the local coordinate system and the global coordinate system is used. Transform the target surface node degrees of freedom to the local contact coordinate system:

[0235] (twenty two)

[0236] in: .

[0237] Assuming the driven wheel transmits the load applied by the driving wheel through tooth meshing, and that it is not subject to external loads, the external load on the target surface is set here. It is zero.

[0238] It is the generalized column vector of the contact forces acting on the target surface. Generally, the target point... It is not a grid node.

[0239] According to the principle of virtual work, the contact force at the four nodes on the target unit is... Along virtual displacement The work done in vain is equal to the target point. Contact force Along virtual displacement The work done in name only. The target point. Contact force Along virtual displacement The work done in the virtual way can be expressed as follows:

[0240] (twenty three)

[0241] Summing all target nodes of the corresponding target element, the total virtual work can be summarized as follows:

[0242] (twenty four)

[0243] Contact forces at four nodes on the target unit Along virtual displacement The work done in the virtual way can be expressed as follows:

[0244] (25)

[0245] By comparing equations (24) and (25), the relationship between nodal forces and non-nodal forces can be obtained:

[0246] (26)

[0247] in: ;

[0248] .

[0249] Substituting equations (21) and (26) back into the equilibrium equation (11) for the potential target surface of the driven wheel, we obtain:

[0250] (27)

[0251] Multiply both sides of the above equation by We can obtain the expression of the target surface equilibrium equation in the local coordinate system because .

[0252] (28)

[0253] in: The other steps and parameters are the same as those in Specific Implementation Methods One through Eight.

[0254] Specific Implementation Method Ten: The difference between this implementation method and Specific Implementation Methods One through Nine is that:

[0255] In S6.6, a second-order conical linear complementary model is constructed by combining the equilibrium equations of the driving wheel and the driven wheel in the local coordinate system; the specific process is as follows:

[0256] S6.6.1: By integrating the equilibrium equations of the driving wheel and the driven wheel in the local coordinate system, we obtain the integrated equilibrium equation in the local coordinate system, which can be expressed as:

[0257] (29)

[0258] Since the frictional contact forces between the two tooth surfaces satisfy Newton's second law, meaning the sum of the contact force vectors at the contact points equals zero, this can be expressed by the formula:

[0259] .

[0260] S6.6.2: Based on the interaction force relationship, the integrated equilibrium equations in the local coordinate system are transformed for the first time. The integrated equilibrium equations in the local coordinate system after the first transformation are expressed by the following formula:

[0261] (30)

[0262] In the formula, For unit array;

[0263] S6.6.3: Because the contact state equation of the tooth surface friction model is about the relationship between relative displacements, the nodal displacement vectors in the integrated equilibrium equations in the transformed local coordinate system are converted into the following form:

[0264] (31)

[0265] A second transformation is performed on the integrated equilibrium equations in the local coordinate system after the first transformation. The integrated equilibrium equations in the local coordinate system after the second transformation are expressed by the following formula:

[0266] (32)

[0267] S6.6.4: The second static condensation operation is performed using the substructure method, transforming the model equilibrium equations into relationships only related to relative displacement. The integrated equilibrium equations in the local coordinate system after the second transformation are then transformed a third time. The integrated equilibrium equations in the local coordinate system after the third transformation are expressed by the following formula:

[0268] (33)

[0269] The above formula is derived Substituting into the second row of the equilibrium equations, we can derive the following equation:

[0270] (34)

[0271] make ;

[0272] ;

[0273] The final integrated equilibrium equation in the local coordinate system under the third transformation is expressed by the following formula:

[0274] (35)

[0275] S6.6.5: Define the relative displacement vector of the tooth surface friction contact model based on theoretical foundations:

[0276] (36)

[0277] In the formula, It is the tangential slip vector of the contact pair element in the two tangential directions of the local coordinate system. It is the normal gap between contact pairs.

[0278] In the initial unloaded state, along Initial gap for direction measurement Its initial gap It can be expressed based on node coordinate data as follows:

[0279] (37)

[0280] in: It is the contact point Node coordinates, The target point The coordinates of the node; This indicates the number of nodes on the potential contact surface of the driving wheel;

[0281] set up This is the generalized relative displacement vector of the two friction contact surfaces after they reach equilibrium at the end of the final load step, including the normal clearance vector and the tangential slip vector.

[0282] In the integrated equilibrium equations defined under the local coordinate system after the third transformation... equal to the initial gap vector and relative displacement vector The difference is expressed by the formula:

[0283] (38)

[0284] In the formula, This represents the generalized relative displacement vector of the two friction contact surfaces after they reach equilibrium at the end of the final load step.

[0285] in,

[0286] in: ;

[0287] ;

[0288] Substituting equation (38) into equation (35), the integrated equilibrium equation in the local coordinate system after the third transformation is transformed for the fourth time. The integrated equilibrium equation in the local coordinate system after the fourth transformation is expressed by the following formula:

[0289] (39)

[0290] S6.6.6: Contact force vector in the integrated equilibrium equations under the local coordinate system of the fourth transformation and relative displacement vector The derivation is in the form of a second-order conic linear complementary equation, expressed by the formula:

[0291]

[0292]

[0293] in, Represents the constraint matrix. Represents the potential matrix, Representing the form of second-order conic linear complementary equations , Representing the form of second-order conic linear complementary equations This can be expressed as a formula:

[0294]

[0295] ,

[0296] In the formula, Represents the k-th element in the potential matrix. Let k represent the k-th element in the constraint matrix; the k-th element in the potential matrix and the k-th element in the constraint matrix are expressed by the following formula:

[0297] (40)

[0298] (41)

[0299] The second-order conic linear complementary equations are in the form of and Substituting the integrated equilibrium equation (39) into the local coordinate system after the fourth transformation, and performing the fifth transformation on the integrated equilibrium equation in the local coordinate system after the fourth transformation, we obtain the formula for the second-order conic linear complementary model as follows:

[0300] (42)

[0301] Finally, a second-order conical linear complementary model of point-to-surface friction contact on the tooth surface is obtained:

[0302]

[0303] In the formula, This represents the first intermediate matrix. The second intermediate matrix is ​​represented by the formula:

[0304]

[0305]

[0306]

[0307]

[0308] (43)

[0309] In the formula, Represents a three-dimensional second-order cone. This represents a one-dimensional non-negative cone. The corresponding vector is... and The problem pertains to second-order cones, and since both components satisfy the complementarity condition, it can be solved using mathematical programming algorithms. The constructed matrix and cone are input into a mathematical computing environment (such as MATLAB). A smooth Newton algorithm, specifically designed for solving complementarity problems, is employed. This algorithm transforms the complementarity problem into a smooth system of equations using a smooth function (such as the Fischer-Burmeister function), and then solves it rapidly using Newton's iteration method, exhibiting superlinear convergence.

[0310] Results Output and Verification: The solver directly outputs the final state of all contact pairs, including: normal contact force, tangential friction force, normal clearance, and tangential slip. These results can be further post-processed to plot tooth surface contact stress contour maps, friction stress distribution, contact imprints, etc., and compared and verified with traditional finite element results.

[0311] The other steps and parameters are the same as those in any of the specific implementation methods one to nine.

[0312] The modeling flowchart for second-order conical linear complementary tooth surface friction contact is as follows: Figure 4 As shown. The detailed process of the second-order conical complementary model for tooth surface friction is as follows:

[0313] (1) Input the gear geometry parameters, use Matlab to numerically discretize the analytical expressions of the involute, addendum circle, dedendum circle and transition curve, and export the discrete node coordinates of the tooth profile curve in TXT text format.

[0314] (2) Use the node connection method to connect the planned tooth profile nodes in sequence to generate tooth profile end face elements. Use the VDRAG stretch command to complete the establishment of the solid element in the tooth width direction of a single tooth. Then use the VGEN rotation complex command on the single tooth element to complete the establishment of the entire gear model.

[0315] (3) According to the gear meshing principle, the initial meshing phase of the two gears is set to complete the assembly. The driven gear is fully constrained, and the driving gear is constrained except for the axial rotational degree of freedom and the working torque is applied.

[0316] (4) The potential contact area of ​​the two gears is refined based on the tooth surface refinement method. This part of the nodes is selected as the main degree of freedom of the model, and the node coordinate data after refinement is exported.

[0317] (5) Using the substructure technology in ANSYS, the model is condensed into a super element (MATRIX50 matrix element) with the potential contact area node as the main node, and the finite element substructure is solved. The stiffness matrix and load array of the contact surface and target surface are exported using the / DEBUG command.

[0318] (6) Use Python to clean the stiffness matrix text and export the stiffness matrix in .mat format. Substitute the coordinate data and stiffness matrix into the second-order conical linear complementary model of tooth surface friction. The complementary equation can be solved by the smooth Newton algorithm based on orthogonal functions or the smooth Newton algorithm based on single-parameter smooth functions.

[0319] (7) Finally, the calculation results of tooth surface normal contact force, normal clearance, tangential friction force and tangential slip are post-processed and displayed.

[0320] The implementation and results of this invention are described in detail below using a pair of parallel-axis herringbone gears as an example. A numerical example analysis is performed based on the modeling process of the above-described tooth surface friction complementary model. Gear geometric parameters and material properties are defined, and a working load is applied at the point where the driving gear's wheel hole is coupled to the center node. The driven gear's wheel hole is coupled to the center node, and all degrees of freedom are constrained. The friction coefficient μ is set to 0.3. The specific gear parameters are shown in Table 1 below.

[0321] Table 1. Basic parameters of the herringbone gear in the example analysis

[0322]

[0323] The solution time of the finite element contact analysis was recorded in ANSYS. Then, the solution time of the second-order conical complementarity method was decomposed into ANSYS substructure stiffness matrix extraction, stiffness matrix and node data processing, and Matlab complementarity program calculation. The working time of each of the three steps was recorded to obtain the total time of the second-order conical complementarity model. The solution time of the two was compared as shown in Table 2. The total solution time of ANSYS software is 4.22 times that of the second-order conical complementarity model, which proves that the second-order conical complementarity model constructed in this invention has a great advantage in solution efficiency.

[0324] Table 2 Comparison of computation time for different methods of herringbone tooth contact analysis

[0325]

[0326] As shown in Table 3, the maximum relative error of the two methods along the entire tooth width (including the stress concentration area) is within 10%. If the stress concentration effect of the tooth tip and tooth root contact is ignored, the relative error of the two methods will be even lower, indicating that the second-order cone model can guarantee the solution accuracy.

[0327] Table 3. Relative errors of the two methods for solving contact stress

[0328]

[0329] like Figure 5 The figure shows the distribution of contact stress projected onto the meshing plane (the directions of the internal common tangents of the two base circles of the gear, i.e., the meshing line direction D and the tooth width direction B) obtained from the second-order cone model. Stress concentration is ignored; the contact stress on the tooth surface is uniformly distributed along the contact line, and each section of the stress surface exhibits a parabolic shape according to Hertzian contact theory. Figure 6 and Figure 7 The figure shows the distribution of frictional stress projected onto the meshing plane by the solution of the second-order cone model. It has the same distribution range as the normal contact stress, and each cross section of the frictional stress surface also presents a parabolic shape.

[0330] (1) Efficiency comparison: The total time consumed in each step of the method of the present invention (data processing, solution, etc.) was recorded and compared with the total time consumed in the standard nonlinear contact analysis of ANSYS. The results are shown in Table 2. The total time consumed by the present invention (95 minutes and 40 seconds) is much lower than that of ANSYS (403 minutes and 59 seconds).

[0331] (2) Accuracy Comparison: The contact stress curves in the tooth width direction calculated by the present invention and ANSYS were extracted and compared. The relative error of the maximum contact stress of several key tooth pairs was calculated, as shown in Table 3. The maximum error was 9.99%, and the average error was within the acceptable range, which fully verified the accuracy of the present invention.

[0332] A statically condensed finite element model of the gear transmission system was constructed, condensing the model into a substructure model with the nodes of the potential contact area on the tooth surface as the main degrees of freedom. Point-to-surface frictional contact elements were defined, and a point-to-surface complementary second-order conical model of tooth surface friction was derived by combining the second-order conical expression of the frictional contact control equation and the principle of virtual work. A numerical example was performed on two parallel-shaft herringbone gears to analyze frictional contact. Comparison with the calculation results from ANSYS software showed that the maximum relative error between the two methods along the entire tooth width (including the stress concentration area) was within 10%. If the stress concentration effect at the tooth tip and root contact is ignored, the relative error would be even lower, demonstrating the accuracy of the second-order conical model. Furthermore, the total solution time of ANSYS software was 4.22 times that of the complementary second-order conical model, demonstrating the high efficiency of the complementary second-order conical model.

[0333] The above description is merely of preferred embodiments of the present invention. It should be understood that the present invention is not limited to the specific embodiments described above. Although the present invention has been disclosed above with reference to preferred embodiments, it is not intended to limit the present invention. Any person skilled in the art can make some modifications or alterations to the above-disclosed technical content to create equivalent embodiments without departing from the scope of the present invention. Any simple modifications, equivalent substitutions, and improvements made to the above embodiments without departing from the scope of the present invention, based on the technical essence of the present invention, and within the spirit and principles of the present invention, shall still fall within the protection scope of the present invention.

Claims

1. A method for calculating tooth surface friction based on mathematical programming algorithms, characterized in that, include: S1: Obtain the geometric parameter data of the driving gear and the driven gear to be calculated; Based on the geometric parameter data of the driving gear and the driven gear to be calculated, mathematical calculation software is used to calculate the tooth profile node coordinate data of the driving gear and the tooth profile node coordinate data of the driven gear. S2: Based on the tooth profile node coordinate data of the driving gear and the tooth profile node coordinate data of the driven gear, parametric modeling is performed using finite element software to obtain the finite element model of the driving gear and the finite element model of the driving gear. S3: In the finite element software, the finite element models of the driving gear and the driving gear are assembled according to the meshing principle to obtain the gear pair meshing finite element model. S4: Determine the potential contact area of ​​the gear pair meshing finite element model, refine the potential contact area using the grid method in the finite element software, and define the nodes of the refined contact area as the master degree of freedom nodes of the tooth surface of the gear pair meshing finite element model. S5: In finite element software, the substructure static condensation method is used to condense all degrees of freedom in the gear pair meshing finite element model onto the main degree of freedom nodes on the tooth surface to obtain a super element model. The super-element model includes: a driving wheel super-element model and a driven wheel super-element model; Export the stiffness matrix of the driving wheel, the first load vector of the driving wheel, the second load vector of the driving wheel, and the nodal coordinate information of the driving wheel in the driving wheel super-element model; Export the driven wheel stiffness matrix, the first load vector of the driven wheel, the second load vector of the driven wheel, and the nodal coordinate information of the driven wheel in the driven wheel super-element model; S6: Based on the stiffness matrix of the driving wheel, the first load vector of the driving wheel, the second load vector of the driving wheel, the nodal coordinate information of the driving wheel, the stiffness matrix of the driven wheel, the first load vector of the driven wheel, the second load vector of the driven wheel, and the nodal coordinate information of the driven wheel, construct a point-plane second-order conical linear complementary model; S7: The second-order conical linear complementary model is solved using the smooth Newton algorithm in mathematical calculation software to obtain the tooth surface friction calculation results.

2. The tooth surface friction calculation method based on mathematical programming algorithm according to claim 1, characterized in that, The gear geometric parameter data to be calculated in S1 includes: involute data, addendum circle data, and dedendum circle data; In step S1, based on the geometric parameter data of the driving gear and the driven gear to be calculated, mathematical calculation software is used to calculate the tooth profile node coordinate data of the driving gear and the driven gear. The specific process is as follows: Input the geometric parameters of the gear to be calculated and the geometric parameters of the driven gear to be calculated into the mathematical calculation software; Based on the geometric parameters of the driving gear to be calculated, the parameterized modeling function of mathematical calculation software is used to discretely generate the tooth profile node coordinate data of the driving gear. Based on the geometric parameters of the driven gear to be calculated, the parameterized modeling function of mathematical calculation software is used to discretely generate the tooth profile node coordinate data of the driven gear.

3. The tooth surface friction calculation method based on mathematical programming algorithm according to claim 2, characterized in that, In step S2, based on the tooth profile node coordinate data of the driving gear and the driven gear, parametric modeling is performed using finite element software to obtain the finite element models of the driving gear and the driven gear; the specific process is as follows: Input the tooth profile node coordinate data of the driving gear and the driven gear obtained from S1 into the finite element software; Based on the tooth profile node coordinate data of the driving gear, the node-element direct generation function in the finite element software is used to generate the finite element model of the driving gear. Based on the tooth profile node coordinate data of the driven gear, the finite element model of the driven gear is generated using the node-element direct generation function in the finite element software.

4. The tooth surface friction calculation method based on mathematical programming algorithm according to claim 3, characterized in that, In step S3, the finite element models of the driving gear and the driven gear are assembled according to the meshing principle to obtain the gear pair meshing finite element model. The specific process is as follows: S3.1: Based on the gear meshing principle, set the initial meshing phase of the driving gear and the initial meshing phase of the driven gear; S3.2: Set the constraints for the finite element model of the driven gear. S3.3: Set the constraints of the finite element model of the driving gear to obtain the meshing finite element model of the gear pair.

5. The tooth surface friction calculation method based on mathematical programming algorithm according to claim 4, characterized in that, In step S6, a point-plane second-order conical linear complementary model is constructed based on the stiffness matrix of the driving wheel, the first load vector of the driving wheel, the second load vector of the driving wheel, the node coordinate information of the driving wheel, the stiffness matrix of the driven wheel, the first load vector of the driven wheel, the second load vector of the driven wheel, and the node coordinate information of the driven wheel; the specific process is as follows: S6.1: Based on the coordinate information of the driving wheel node and the driven wheel node, K point-to-surface friction contact pairs are constructed in the mathematical calculation software. Each point-to-surface friction contact pair consists of a contact node and a target surface. The contact node is a node on the active gear tooth surface, and the target surface is a four-node isoparametric element on the driven gear tooth surface; S6.2: Establish a local coordinate system for each point-to-surface friction contact pair element to obtain K local coordinate systems; And define a generalized contact force vector and a generalized relative displacement vector in each local coordinate system; S6.3: Construct the transformation matrix based on the K local coordinate systems; S6.4: Based on the transformation matrix, the stiffness matrix of the driving wheel, the first load vector of the driving wheel, and the second load vector of the driving wheel, construct the equilibrium equation of the driving wheel in the local coordinate system; S6.5: Based on the transformation matrix, the driven wheel stiffness matrix, the driven wheel first load vector, and the driven wheel second load vector, construct the driven wheel equilibrium equation in the local coordinate system; S6.6: Construct a second-order conical linear complementary model by combining the equilibrium equations of the driving wheel in the local coordinate system and the equilibrium equations of the driven wheel in the local coordinate system.

6. The tooth surface friction calculation method based on mathematical programming algorithm according to claim 5, characterized in that, In step S6.2, a local coordinate system is established for each point-to-surface friction contact pair, resulting in K local coordinate systems; the specific process is as follows: To establish a local coordinate system for the k-th point-to-surface friction contact pair element, the specific process for obtaining the k-th local coordinate system is as follows: Let the contact nodes in the k-th point-to-surface friction contact element be at a certain direction angle. Along Move in the direction normal, and in the direction relative to the target surface. and Tangential slip, , and Both are unit vectors. According to unit vector Unit vector With unit vector This forms the k-th local spatial coordinate system.

7. The tooth surface friction calculation method based on mathematical programming algorithm according to claim 6, characterized in that, In S6.3, the transformation matrix is ​​constructed based on K local coordinate systems, which can be expressed by the following formula: In the formula, Represents the transformation matrix. express The element in the first row and first column, express The element in the k-th row and k-th column, express The element in the Kth row and Kth column, the The element in row k and column k The calculation formula is expressed as: In the formula, x represents the x-direction vector, y represents the y-direction vector, and z represents the z-direction vector. This represents the cosine value.

8. The tooth surface friction calculation method based on mathematical programming algorithm according to claim 7, characterized in that, In step S6.4, the equilibrium equations of the driving wheel in the local coordinate system are constructed based on the transformation matrix, the stiffness matrix of the driving wheel, the first load vector of the driving wheel, and the second load vector of the driving wheel; expressed by the formula: In the formula, This represents the stiffness matrix of the driving wheel in the local coordinate system. This represents the first load vector of the driving wheel in the local coordinate system. This represents the second load vector of the driving wheel in the local coordinate system. The displacement vector of the driving wheel node in the local coordinate system is expressed by the formula: In the formula, express The inverse matrix, Represents the stiffness matrix of the drive wheel, Indicates the first load vector of the driving wheel, This represents the second load vector of the driving wheel; This represents the displacement vector of the active wheel node.

9. The tooth surface friction calculation method based on mathematical programming algorithm according to claim 8, characterized in that, In S6.5, based on the transformation matrix, the driven wheel stiffness matrix, the driven wheel first load vector, and the driven wheel second load vector, the driven wheel equilibrium equation in the local coordinate system is constructed, which is expressed by the following formula: In the formula, This represents the stiffness matrix of the driven wheel in the local coordinate system. This represents the displacement vector of the driven wheel node in the local coordinate system. The second load vector of the driven wheel in the local coordinate system is expressed by the formula: In the formula, This represents the stiffness matrix of the driven wheel. This represents the displacement vector of the driven wheel node. This represents the second load vector of the driven wheel. Represents a pointer matrix. Represents a shape function matrix, pointer matrix Expressed as a formula: In the formula, express The first element, express The kth element, express The Kth element, Shape function matrix Expressed as a formula: express, express The element in the first row and first column, express The element in the k-th row and k-th column, express The element in row K and column K.

10. The tooth surface friction calculation method based on mathematical programming algorithm according to claim 9, characterized in that, In S6.6, a second-order conical linear complementary model is constructed by combining the equilibrium equations of the driving wheel and the driven wheel in the local coordinate system; the specific process is as follows: S6.6.1: By integrating the equilibrium equations of the driving wheel and the driven wheel in the local coordinate system, we obtain the integrated equilibrium equation in the local coordinate system, which can be expressed as: （29） S6.6.2: Perform the first transformation on the integrated equilibrium equations in the local coordinate system. The integrated equilibrium equations in the local coordinate system after the first transformation are expressed by the following formula: （30） In the formula, For unit array; S6.6.3: Perform a second transformation on the integrated equilibrium equations in the local coordinate system after the first transformation. The integrated equilibrium equations in the local coordinate system after the second transformation are expressed by the following formula: （32） S6.6.4: The integrated equilibrium equations in the local coordinate system after the second transformation are subjected to a third transformation. The integrated equilibrium equations in the local coordinate system after the third transformation are expressed by the following formula: （35） S6.6.5: In the integrated equilibrium equations under the local coordinate system of the third transformation, equal to the initial gap vector and relative displacement vector The difference is expressed by the formula: （38） In the formula, This represents the generalized relative displacement vector of the two friction contact surfaces after they reach equilibrium at the end of the final load step. in, The integrated equilibrium equations in the local coordinate system after the third transformation are subjected to a fourth transformation. The integrated equilibrium equations in the local coordinate system after the fourth transformation are expressed by the following formula: （39） S6.6.6: In the integrated equilibrium equations under the local coordinate system of the fourth transformation and The derivation is in the form of a second-order conic linear complementary equation, expressed by the formula: in, Represents the constraint matrix. Represents the potential matrix, Representing the form of second-order conic linear complementary equations , Representing the form of second-order conic linear complementary equations This can be expressed as a formula: ， In the formula, Represents the k-th element in the potential matrix. This represents the k-th element in the constraint matrix; A fifth transformation is performed on the integrated equilibrium equations in the local coordinate system after the fourth transformation, yielding the formula for the second-order conic linear complementary model: In the formula, This represents the first intermediate matrix. The second intermediate matrix is ​​represented by the formula: 。

Images