Gear pair meshing stiffness calculation method considering thermoelastic coupling effect

Abstract

The invention discloses a gear pair meshing stiffness calculation method considering a thermoelastic coupling effect, which relates to the field of gear transmission, and comprises the following steps: establishing a gear temperature field model and a contact model by utilizing a finite element method, obtaining a gear heat conduction matrix and a structural stiffness matrix, and calculating to obtain gear meshing stiffness and tooth surface load distribution based on an actual tooth surface equation. Calculating tooth surface heat flux density based on tooth surface load distribution, tooth surface friction coefficient and relative sliding speed of each meshing position to obtain a gear thermal load matrix, calculating to obtain an equivalent node load matrix according to a gear heat conduction matrix and the thermal load matrix, and calculating to obtain tooth surface thermal deformation by combining the equivalent node load matrix with a gear structure stiffness matrix. And repeating the steps to obtain the meshing rigidity of the gear pair under the thermoelastic coupling effect. According to the method, the gear contact state iterative solution model based on the potential energy method, the matrix polycondensation method and the slicing thought is established, the requirement for computing equipment is lowered, and the solution efficiency is improved.

Description

Technical Field

[0001] This invention relates to the field of gear transmission, and in particular to a method for calculating the meshing stiffness of gear pairs considering the thermo-elastic coupling effect. Background Technology

[0002] Gear drives are widely used in aerospace, wind power generation, and shipbuilding, and their contact state and dynamic characteristics directly affect the reliability and safety of mechanical devices. As an important internal excitation source, gear meshing stiffness analysis is crucial for predicting and optimizing system vibration and noise. For high-power-density gear drive systems, the transmission power is large, the operating speed is high, and the system temperature rise is large. There is a strong coupling effect between the gear contact state and thermal characteristics, which in turn affects the gear meshing stiffness. However, current gear thermo-elastic coupling analysis is usually carried out using commercial software, which has a high computational cost. Therefore, it is necessary to propose a rapid analysis method for gear pair meshing stiffness that considers thermo-elastic coupling.

[0003] To address these issues, a method for calculating the meshing stiffness of gear pairs that considers the thermo-elastic coupling effect is urgently needed. Summary of the Invention

[0004] To address the aforementioned issues, this application proposes a method for calculating the meshing stiffness of gear pairs that considers the thermoelastic coupling effect, thereby resolving the problems existing in the prior art.

[0005] A method for calculating the meshing stiffness of a gear pair considering the thermo-elastic coupling effect includes the following steps:

[0006] S1. Establish the gear temperature field model and contact model using the finite element method (ANSYS software), and extract the gear heat conduction matrix [K] from the gear temperature field model and contact model. T and structural stiffness matrix [K] S ];

[0007] S2. Based on the actual tooth surface equation, the gear meshing stiffness K(t) and tooth surface load distribution F are calculated. mn (t);

[0008] S3, Obtain the tooth surface friction coefficient μ mn (t), relative sliding velocity v at each meshing position mn (t), based on the tooth surface load distribution, tooth surface friction coefficient, and relative sliding velocity at each meshing position, the tooth surface heat flux density is calculated, and the gear heat load matrix [F] is obtained. T ];

[0009] S4. According to the gear heat conduction matrix [K] T ] and thermal load matrix [F T ], calculate the gear node temperature matrix [T q], and thus the equivalent nodal load matrix [F] is obtained. TS ];

[0010] Combining the equivalent nodal load matrix [F TS ] and gear structure stiffness matrix [K S The thermal deformation of the tooth surface [x] was calculated. q (T q )]、[y q (T q )];

[0011] S5. Repeat steps S2 to S4 until the relative error of the gear pair meshing stiffness meets the iterative convergence condition, and obtain the gear pair meshing stiffness under thermo-elastic coupling.

[0012] Preferably, the gear structure stiffness matrix [K] is extracted. S Then, the Guyan matrix reduction method is used to obtain the load-displacement relationship matrix of the tooth surface nodes, which reduces the calculation of tooth surface thermal deformation and then calculates the amount of tooth surface thermal deformation.

[0013] The load-displacement relationship matrix of the original nodes is as follows:

[0014]

[0015] The load-displacement relationship matrix of the tooth surface nodes after reduction using the Guyan matrix can be expressed as:

[0016]

[0017] In the formula, a represents the number of the reducible node set, b represents the number of the tooth surface node set, and K aa To reduce the principal stiffness matrix region of the nodes, K ab For the right-hand matrix region corresponding to the shrinkable principal stiffness matrix region of the nodes, K ba To reduce the lower matrix region corresponding to the principal stiffness matrix region of the nodes, K bb X represents the stiffness matrix region corresponding to the tooth surface nodes. a X is the displacement matrix corresponding to the non-tooth surface nodes. b F is the displacement matrix corresponding to the tooth surface nodes. a F is the load matrix corresponding to the non-tooth surface nodes. b This is the load matrix corresponding to the tooth surface nodes. The equivalent stiffness matrix after the reduction of the Guyan matrix. Let I be the equivalent load matrix after the reduction of the Guyan matrix, and let I be the identity matrix. For K aa The inverse matrix of , where T is the matrix transpose symbol.

[0018] Preferably, the actual tooth surface equation in S2 is based on the machining principle of the hypothetical machining tool. The tooth profile is modified by replacing the original tooth profile curve at the tooth tip and tooth root positions of the hypothetical machining tool end face. The tooth profile is modified by controlling the axial movement of the grinding wheel. The tooth profile is modified by combining the tooth surface thermal deformation.

[0019] Preferably, the expression for the actual tooth surface equation is:

[0020]

[0021] Among them, l z α is the distance from any point on the tooth surface to the gear end face; c1 is the major semi-axis of the ellipse arc, d1 is the minor semi-axis of the ellipse arc; θ is the angular parameter at the tip of the ellipse; α t β is the gear end face pressure angle, and β is the gear helix angle; The rotation angle is r; the pitch circle radius of the gear being machined is r. Here is the equation of the tooth surface in coordinate system S3. Let x be the x-direction component of the tooth surface equation in coordinate system S2. Let C be the y-direction component of the tooth surface equation in coordinate system S2. C x represents the tooth profile modification at any point on the gear tooth. q Let y be the thermal deformation of the gear teeth in the x direction. q T represents the thermal deformation of the gear teeth in the y-direction. q R represents the gear pitch temperature, i represents the number of tooth surface segments, and R represents the gear pitch temperature. S1 Let x be the equation of the tooth surface in coordinate system S1. D Let x be the x-coordinate of point D in coordinate system S1. C Let x be the x-coordinate of point C in coordinate system S1. Oe Let x be the x-coordinate of the ellipse center in coordinate system S1, u1 be the tooth tip modification amount, u2 be the tooth root modification amount, a1 be the tooth tip modification curve coefficient, a2 be the tooth root modification curve coefficient, n be the order of the modification curve, and y be the x-coordinate of the ellipse center in coordinate system S1. D Let y be the y-coordinate of point D in coordinate system S1. C Let y be the y-coordinate of point C in coordinate system S1.

[0022] Preferably, in S2, the potential energy method and the slicing concept are used to calculate the gear meshing stiffness K(t) and the tooth surface load distribution F. mn (t);

[0023] Based on the slicing concept, the gear is sliced ​​along the tooth width direction. When calculating the time-varying meshing stiffness, each slice of tooth is simplified into a cantilever beam model. The end face deformation, axial deformation and base deformation of the tooth are comprehensively considered to obtain the meshing stiffness of the gear pair.

[0024] Based on the actual tooth surface equation, solve for the tooth profile deviation δ of each thin-plate gear tooth at time t.mn (t) and meshing stiffness k mn (t).

[0025] Preferably, if the sum of the loads borne by each thin-plate tooth at time t is equal to the load F transmitted by the gear pair, then the following should be satisfied:

[0026]

[0027] F mn (t)=(δ(t)-δ mn (t))k mn (t);

[0028] In the formula, m is the number of gear pairs participating in meshing at time t; n is the number of thin-plate teeth on each meshing gear pair; t is time; F is the load transmitted by the gear pair; K(t) is the meshing stiffness of the gear pair at time t; δ(t) is the deformation of the gear teeth under load at time t; δ mn (t) represents the profile deviation of each thin-plate gear tooth at time t.

[0029] Preferably, the expression for meshing stiffness is:

[0030]

[0031] In the formula: j = 1, 2 represent the driving and driven gears; k tb k is the end face bending stiffness. ts For end face shear stiffness, k ta For radial compressive stiffness, k tf k is the radial stiffness of the matrix. ab For axial bending stiffness, k at For axial torsional stiffness, k af k is the axial stiffness of the matrix. h Here, is the Hertzian contact stiffness, and k is the meshing stiffness.

[0032] Preferably, the tooth surface friction coefficient μ in S3 mn (t) is calculated using load sharing theory and elastohydrodynamic lubrication theory;

[0033] The relative sliding velocity at the engagement position v mn (t) is calculated using the tooth surface slip ratio formula. Preferably, the expression for the comprehensive friction coefficient of the tooth surface is:

[0034] μ mn (t)=γ(t)μ a (t)+(1-γ(t))μ b ;

[0035] In the formula, μ a (t) is the oil film friction coefficient, calculated based on elastohydrodynamic lubrication theory; μb The roughness peak friction coefficient is determined based on material properties and roughness; t is time, μ mn γ(t) is the comprehensive friction coefficient at time t, and γ(t) is the oil film bearing ratio at time t.

[0036] Preferably, the calculation method for the oil film bearing ratio includes:

[0037] Define the initial oil film load ratio γ(t=0=1); perform iterative solutions using a thermo-elasto-hydrodynamic lubrication analysis process; until the rough peak contact pressure p is calculated according to the following formula. c1 and p c2 If the difference is less than the allowable value, the actual oil film bearing ratio is obtained.

[0038]

[0039] In the formula, E0 is the equivalent elastic modulus of the contact body, R is the combined radius of curvature of the contact point, L is the contact line length of the tooth surface, and β s Where n is the radius of the rough peak. d For the rough peak density, σ s d d The roughness is a highly correlated value, where d d =1.15σ s h c Where W is the thickness of the oil film at the center, and F is the dimensionless load. 5 / 2 This is a function corresponding to changes in film thickness.

[0040] Preferably, in S4, the expression for the amount of thermal deformation of the tooth surface is:

[0041]

[0042] [T q ] = [K T ] -1 [F T ];

[0043] Where, x q Let y be the thermal deformation of the gear teeth in the x direction. q Let z be the thermal deformation of the gear teeth in the y-direction. q T represents the thermal deformation of the gear teeth in the z-direction. q For gear pitch temperature, K S Let F be the stiffness matrix of the gear structure. TS For the equivalent nodal load matrix, [K S ] -1 For K S The inverse matrix.

[0044] In summary, the gear pair meshing stiffness calculation method considering the thermo-elastic coupling effect of the present invention has the following advantages compared with traditional technologies:

[0045] 1. Taking into account tooth profile modification, thermal deformation, and operating conditions, the meshing stiffness of gear pairs considering the thermo-elastic bidirectional coupling effect can be calculated quickly.

[0046] 2. The load-displacement relationship matrix of the tooth surface nodes can be obtained by using the stiffness reduction method, which can reduce the calculation cost of gear thermo-elastic bidirectional coupling analysis and improve the solution efficiency.

[0047] 3. There is no need to establish a complex gear thermo-elastic two-way coupling simulation model. By establishing an iterative solution model for gear contact state based on the potential energy method, matrix contraction method and slicing idea, the requirements for computing equipment are reduced and the solution efficiency is improved.

[0048] The technical method of the present invention will be further described in detail below with reference to the accompanying drawings and embodiments. Attached Figure Description

[0049] Figure 1 This is a flowchart of the steps in the method for calculating the meshing stiffness of a gear pair considering the thermo-elastic coupling effect of the present invention.

[0050] Figure 2 This invention provides a method for calculating the meshing stiffness of a gear pair considering thermo-elastic coupling effects, which includes a gear solid model and a discretized mesh model.

[0051] Figure 3 This is a flowchart of the solution for gear meshing stiffness and tooth surface load in a gear pair meshing stiffness calculation method considering thermo-elastic coupling effect according to the present invention.

[0052] Figure 4 This is a schematic diagram illustrating the derivation of the tooth surface equation in a method for calculating the meshing stiffness of a gear pair considering the thermo-elastic coupling effect, as described in this invention. Figure 4 (a) in the diagram is a schematic diagram of the hypothetical tool tooth profile. Figure 4 (b) in the diagram is a schematic diagram of the kinematic relationship between the hypothetical cutting tool and the gear being machined. Figure 4 (c) in the diagram is a schematic diagram of the kinematic relationship between the grinding wheel and the gear being machined. Figure 4 (d) in the diagram is a schematic diagram of the axial feed motion of the grinding wheel;

[0053] Figure 5 This is the result of gear contact analysis calculation in the gear pair meshing stiffness calculation method considering thermo-elastic coupling effect of the present invention. Figure 5 In the diagram, (a) represents the load distribution on the tooth surface. Figure 5 In the figure, (b) represents the overall friction coefficient of the tooth surface. Figure 5 In this context, (c) represents the relative sliding velocity of the tooth surfaces;

[0054] Figure 6 The gear meshing stiffness calculation result is shown in the gear pair meshing stiffness calculation method considering the thermo-elastic coupling effect of the present invention. Detailed Implementation

[0055] The technical method of the present invention will be further described below with reference to the accompanying drawings and embodiments. It should be noted that, unless otherwise specifically stated, the relative arrangement, numerical expressions, and values ​​of the components and steps described in these embodiments do not limit the scope of this application.

[0056] The following description of at least one exemplary embodiment is merely illustrative and is in no way intended to limit the scope of this application and its application or use.

[0057] Techniques, systems, and equipment known to those skilled in the art may not be discussed in detail, but where appropriate, they should be considered part of the instruction manual.

[0058] In all the examples shown and discussed herein, any specific values ​​should be interpreted as merely exemplary and not as limitations. Therefore, other examples of exemplary embodiments may have different values.

[0059] Unless otherwise defined, the technical or scientific terms used in this invention shall have the ordinary meaning as understood by one of ordinary skill in the art to which this invention pertains.

[0060] This invention provides a method for calculating the meshing stiffness of gear pairs considering the thermo-elastic coupling effect, such as... Figure 1 As shown, it includes the following steps:

[0061] S1. Using ANSYS software, the gear model is discretized into a mesh, such as... Figure 2 As shown, a gear temperature field model and a contact model are established, and the gear heat conduction matrix [K] is extracted. T and structural stiffness matrix [K] S ];

[0062] Furthermore, the gear structure stiffness matrix [K] is extracted. S Then, the Guyan matrix reduction method was used to obtain the load-displacement relationship matrix of the tooth surface nodes, reducing the amount of calculation of tooth surface thermal deformation.

[0063] The load-displacement relationship matrix of the original nodes is as follows:

[0064]

[0065] The load-displacement relationship matrix of the tooth surface nodes after reduction using the Guyan matrix can be expressed as:

[0066]

[0067] In the formula, a represents the number of the reducible node set, b represents the number of the tooth surface node set, and K aa To reduce the principal stiffness matrix region of the nodes, K ab For the right-hand matrix region corresponding to the shrinkable principal stiffness matrix region of the nodes, K ba To reduce the lower matrix region corresponding to the principal stiffness matrix region of the nodes, K bb X represents the stiffness matrix region corresponding to the tooth surface nodes. a X is the displacement matrix corresponding to the non-tooth surface nodes. b F is the displacement matrix corresponding to the tooth surface nodes. a F is the load matrix corresponding to the non-tooth surface nodes. b This is the load matrix corresponding to the tooth surface nodes. The equivalent stiffness matrix after the reduction of the Guyan matrix. Let I be the equivalent load matrix after the reduction of the Guyan matrix, and let I be the identity matrix. For K aa The inverse matrix of , where T is the matrix transpose symbol.

[0068] S2. Based on the actual tooth surface equation, the gear meshing stiffness K(t) and tooth surface load distribution F are obtained by iterative calculation. mn (t), Figure 3 The iterative solution process is presented as follows: Based on the basic parameters, modification parameters, and thermal deformation of the gears, the tooth surface equation is derived. The kinematic relationship and relative positional relationship between the driving and driven gears are determined by combining the installation conditions. Then, the distance between the meshing tooth profiles and the meshing stiffness of each plate gear pair are calculated based on the direction of the gear pair's meshing line. Assuming the meshing deformation at time t, and combining the distance between the meshing tooth profiles, the plate gear pairs participating in the meshing are determined. Based on the meshing stiffness of each plate gear pair, the total load borne by the participating plate gear pairs is calculated by superposition and compared with the load transmitted by the gear pair. If the error between the two is less than the allowable error value, convergence is determined, and the result is output. If the error is greater than the allowable error value, the meshing deformation at time t is modified, and the above operation is repeated until the error between the total load borne by the participating plate gear pairs and the load transmitted by the gear pair is less than the allowable error value.

[0069] Furthermore, considering tooth profile modification and thermal deformation, the actual tooth surface equation is derived. The actual tooth surface equation described in S2 is based on the machining principle of a hypothetical machining tool. Tooth profile modification is achieved by replacing the original tooth profile curve at the tooth tip and root positions on the hypothetical machining tool end face with a tooth profile modification curve. Tooth profile modification is achieved by controlling the axial movement of the grinding wheel. Figure 4 As shown, this is obtained by combining the thermal deformation of the tooth surface.

[0070] Furthermore, the expression for the actual tooth surface equation is as follows:

[0071]

[0072] Among them, l z α is the distance from any point on the tooth surface to the gear end face; c1 is the major semi-axis of the ellipse arc, d1 is the minor semi-axis of the ellipse arc; θ is the angular parameter at the tip of the ellipse; α t β is the gear end face pressure angle, and β is the gear helix angle; The rotation angle is r; the pitch circle radius of the gear being machined is r. Here is the equation of the tooth surface in coordinate system S3. Let x be the x-direction component of the tooth surface equation in coordinate system S2. Let C be the y-direction component of the tooth surface equation in coordinate system S2. C x represents the tooth profile modification at any point on the gear tooth. q Let y be the thermal deformation of the gear teeth in the x direction. q T represents the thermal deformation of the gear teeth in the y-direction. q R represents the gear pitch temperature, i represents the number of tooth surface segments, and R represents the gear pitch temperature. S1 Let x be the equation of the tooth surface in coordinate system S1. D Let x be the x-coordinate of point D in coordinate system S1. C Let x be the x-coordinate of point C in coordinate system S1. Oe Let x be the x-coordinate of the ellipse center in coordinate system S1, u1 be the tooth tip modification amount, u2 be the tooth root modification amount, a1 be the tooth tip modification curve coefficient, a2 be the tooth root modification curve coefficient, n be the order of the modification curve, and y be the x-coordinate of the ellipse center in coordinate system S1. D Let y be the y-coordinate of point D in coordinate system S1. C Let y be the y-coordinate of point C in coordinate system S1.

[0073] Furthermore, in S2, the potential energy method and the slicing concept are used to calculate the gear meshing stiffness K(t) and the tooth surface load distribution F. mn (t);

[0074] Based on the slicing concept, the gear is sliced ​​along the tooth width direction. When calculating the time-varying meshing stiffness, each slice of tooth is simplified into a cantilever beam model. The end face deformation, axial deformation and base deformation of the tooth are comprehensively considered to obtain the meshing stiffness of the gear pair.

[0075] Furthermore, the expression for meshing stiffness is:

[0076]

[0077] In the formula, j = 1 and 2 represent the driving and driven wheels; k tb k is the end face bending stiffness. ts For end face shear stiffness, k ta For radial compressive stiffness, k tf k is the radial stiffness of the matrix.ab For axial bending stiffness, k at For axial torsional stiffness, k af k is the axial stiffness of the matrix. h Here, is the Hertzian contact stiffness, and k is the meshing stiffness.

[0078]

[0079]

[0080] In the formula: x c y c x and y are the coordinates of the meshing points of the teeth of each thin-plate gear; x and y are the coordinates of any point on the gear teeth; x M α is the x-coordinate value at the root of the gear tooth; c β is the end-face pressure angle at the meshing point of each thin-plate gear tooth; β is the gear helix angle; E, G, and μ represent the elastic modulus, shear modulus, and Poisson's ratio of the gear material, respectively; N s L*, M*, P*, and Q* are the number of sections the gear is divided into along the tooth width direction; L*, M*, P*, and Q* are size factors; I af R is the moment of inertia of the cross section at a distance x from the center of the gear; B is the tooth width; r f The radius of the gear tooth root circle is given.

[0081] Pressure angle α at the meshing point of the gear teeth c and the moment of inertia of the cross section I af These can be expressed as:

[0082]

[0083] In the formula: and y is the x and y components of the tooth surface normal at the meshing point of the thin-plate gear teeth, which can be calculated based on the tooth surface equation derived above; r0 is the radius of the gear shaft hole.

[0084] Based on the actual tooth surface equation, solve for the tooth profile deviation δ of each thin-plate gear tooth at time t. mn (t) and meshing stiffness k mn (t). At time t, the sum of the loads borne by each thin-plate tooth is equal to the load F transmitted by the gear pair, then the following should be satisfied:

[0085]

[0086] F mn (t)=(δ(t)-δ mn (t))k mn (t);

[0087] In the formula, m is the number of gear pairs participating in meshing at time t; n is the number of thin-plate teeth on each meshing gear pair; t is time; F is the load transmitted by the gear pair; K(t) is the meshing stiffness of the gear pair at time t; δ(t) is the deformation of the gear teeth under load at time t; δ mn (t) represents the profile deviation of each thin-plate gear tooth at time t.

[0088] S3, Obtain the tooth surface friction coefficient μ mn (t), relative sliding velocity v at each meshing position mn (t), based on tooth surface load distribution, tooth surface friction coefficient, and relative sliding velocity at each meshing position, such as Figure 5 As shown, the heat flux density on the tooth surface is calculated, and the gear thermal load matrix [F] is obtained. T ];

[0089] According to the gear heat conduction matrix [K] T ] and thermal load matrix [F T ], calculate the gear node temperature matrix [T q ], and thus the equivalent nodal load matrix [F] caused by temperature rise is obtained. TS ];

[0090] Furthermore, the tooth surface friction coefficient μ mentioned in S3 mn (t) is calculated using load sharing theory and elastohydrodynamic lubrication theory;

[0091] The relative sliding velocity at the engagement position v mn (t) is calculated using the tooth surface slip ratio formula. Furthermore, the expression for the overall tooth surface friction coefficient is:

[0092] μ mn (t)=γ(t)μ a (t)+(1-γ(t))μ b ;

[0093] In the formula, μ a (t) is the oil film friction coefficient, calculated based on elastohydrodynamic lubrication theory; μ b The roughness peak friction coefficient is determined based on material properties and roughness; t is time, μ mn γ(t) is the comprehensive friction coefficient at time t, and γ(t) is the oil film bearing ratio at time t.

[0094] Furthermore, the calculation methods for the oil film load-bearing ratio include:

[0095] Define the initial oil film load ratio γ(t=0=1); iteratively solve the problem using a thermo-elasto-hydrodynamic lubrication analysis process until the rough peak contact pressure p is calculated according to the following formula. c1 and p c2 If the difference is less than the allowable value, the actual oil film bearing ratio is obtained.

[0096]

[0097] In the formula, E0 is the equivalent elastic modulus of the contact body, R is the combined radius of curvature of the contact point, L is the contact line length of the tooth surface, and β s Where n is the radius of the rough peak. d For the rough peak density, σ s d d The roughness is a highly correlated value, where d d =1.15σ s h c Where W is the thickness of the oil film at the center, and F is the dimensionless load. 5 / 2 This is a function corresponding to changes in film thickness.

[0098] S4. According to the gear heat conduction matrix [K] T ] and thermal load matrix [F T ], calculate the gear node temperature matrix [T q ], and thus the equivalent nodal load matrix [F] is obtained. TS ];

[0099] Equivalent nodal load matrix [F TS Combining the gear structure stiffness matrix [K] S The thermal deformation of the tooth surface [x] was calculated. q (T q )]、[y q (T q )];

[0100] In S4, the expression for the thermal deformation of the tooth surface is:

[0101]

[0102] [T q ] = [K T ] -1 [F T ];

[0103] Where, x q Let y be the thermal deformation of the gear teeth in the x direction. q Let z be the thermal deformation of the gear teeth in the y-direction. q T represents the thermal deformation of the gear teeth in the z-direction. q For gear pitch temperature, K S Let F be the stiffness matrix of the gear structure. TS For the equivalent nodal load matrix, [K S ] -1 For K S The inverse matrix.

[0104] S5. Repeat steps S2 to S4 until the relative error of the gear pair meshing stiffness satisfies the iterative convergence condition, thus obtaining the gear pair meshing stiffness under thermo-elastic coupling. Figure 6 As shown.

[0105] Finally, it should be noted that the above embodiments are only used to illustrate the technical methods of the present invention and not to limit them. Although the present invention has been described in detail with reference to preferred embodiments, those skilled in the art should understand that modifications or equivalent substitutions can still be made to the technical methods of the present invention, and these modifications or equivalent substitutions cannot cause the modified technical methods to deviate from the spirit and scope of the technical methods of the present invention.

Claims

1. A method for calculating mesh stiffness of a gear pair considering thermal-elastic coupling effect, characterized in that, The method comprises the following steps: S1, the finite element method is used to establish gear temperature field model and contact model, extract gear heat conduction matrix [K T ] and structure stiffness matrix [K S ]; S2, based on the actual tooth surface equation, the gear mesh stiffness K(t) and the tooth surface load distribution F are calculated mn (t); S3, acquiring the friction coefficient μ of the tooth surface mn (t), the relative sliding speed v of each meshing position mn (t), based on the tooth surface load distribution, the friction coefficient μ of the tooth surface, and the relative sliding speed v of each meshing position, calculating the tooth surface heat flux density, and obtaining the gear thermal load matrix [F T ] According to the gear heat conduction matrix [K T ] and the thermal load matrix [F T ], the gear node temperature matrix [T q ] is calculated, and then the equivalent node load matrix [F TS ] is obtained. S4, equivalent node load matrix [F TS ] combined with gear structure stiffness matrix [K S ], the tooth surface thermal deformation amount [x q (T q )], [y q (T q )] is calculated; S5, repeating steps S2-S4 until the relative error of the gear pair meshing stiffness meets the iterative convergence condition, and obtaining the gear pair meshing stiffness under the thermal elastic coupling effect.

2. The method for calculating the meshing stiffness of a gear pair considering thermal-elastic coupling effect according to claim 1, characterized in that, Extracting gear structure stiffness matrix [K S ]After, using Guyan matrix condensation method to obtain the tooth surface node load-displacement relationship matrix, and then calculate the tooth surface thermal deformation. The load-displacement relationship matrix of the original node is as follows: The load-displacement relationship matrix of the tooth surface node after adopting the Guyan matrix condensation can be expressed as: where a denotes the reducible node set number, b denotes the tooth surface node set number, K aa is the reducible node main stiffness matrix area, K ab is the right side matrix area corresponding to the reducible node main stiffness matrix area, K ba is the lower side matrix area corresponding to the reducible node main stiffness matrix area, K bb is the tooth surface node corresponding stiffness matrix area, X a is the non-tooth surface node corresponding displacement matrix, X b is the tooth surface node corresponding displacement matrix, F a is the non-tooth surface node corresponding load matrix, F b is the tooth surface node corresponding load matrix, is the Guyan matrix condensed equivalent stiffness matrix, is the Guyan matrix condensed equivalent load matrix, I is the unit matrix, is the inverse matrix of K aa , T is the matrix transpose symbol.

3. The method for calculating the meshing stiffness of a gear pair considering the thermo-elastic coupling effect according to claim 2, characterized in that, The actual tooth surface equation in S2 is based on the machining principle of the imaginary machining tool, adopts a tooth profile modification curve to replace the original tooth profile curve at the addendum and dedendum positions of the end surface tooth profile of the imaginary machining tool to realize tooth profile modification, controls the axial movement of the grinding wheel to realize tooth direction modification, and is combined with tooth surface thermal deformation.

4. The method for calculating the meshing stiffness of a gear pair considering thermal-elastic coupling effect according to claim 3, characterized in that, The expression of the actual tooth surface equation is as follows: wherein, l z is the distance from any point on the tooth surface to the gear end surface; c1 is the long semi-axis of the elliptical arc, d1 is the short semi-axis of the elliptical arc; θ is the angle parameter variable at the tool tip elliptical position; α t is the pressure angle of the gear end surface, and β is the helix angle of the gear; is the rotation angle; r is the pitch circle radius of the machined gear, is the tooth surface equation in the coordinate system S3, is the x-direction component of the tooth surface equation in the coordinate system S2, is the y-direction component of the tooth surface equation in the coordinate system S2, C is the tooth direction modification amount at any point of the gear tooth, x q is the x-direction thermal deformation amount of the gear tooth, y q is the y-direction thermal deformation amount of the gear tooth, T q is the gear node temperature, i is the number of tooth surface segments, R S1 is the tooth surface equation in the coordinate system S1, x D is the x-coordinate value of point D in the coordinate system S1, x C is the x-coordinate value of point C in the coordinate system S1, x Oe is the x-coordinate value of the elliptical center in the coordinate system S1, u1 is the gear tooth tip modification amount, u2 is the gear tooth root modification amount, a1 is the tooth tip modification curve coefficient, a2 is the tooth root modification curve coefficient, and n is the modification curve order, y D is the y-coordinate value of point D in the coordinate system S1, y C is the y-coordinate value of point C in the coordinate system S1.

5. The method for calculating the meshing stiffness of a gear pair considering thermal-elastic coupling effect according to claim 4, characterized in that, In S2, the potential energy method and the slice idea are adopted to calculate the gear mesh stiffness K(t) and the tooth surface load distribution F mn (t); Based on the slicing idea, the gear is sliced along the tooth width direction, each thin slice gear tooth is simplified into a cantilever beam model when the time-varying meshing stiffness is calculated, the end surface deformation, axial deformation and substrate deformation of the gear tooth are comprehensively considered, and the gear pair meshing stiffness is obtained. According to the actual tooth surface equation, the tooth profile deviation δ of each lamella gear tooth at time t is solved mn (t) and the meshing stiffness k mn (t); the sum of the load borne by each lamella gear tooth at time t is equal to the load F transmitted by the gear pair, so that: F mn (t) = (δ(t) - δ mn (t))k mn (t); In the formula, m is the number of tooth pairs of the gear pair participating in engagement at time t; n is the number of sheet gears bearing on each meshing tooth pair, t is time, F is the load transmitted by the gear pair, K(t) is the engagement stiffness of the gear pair at time t, δ(t) is the deformation of the gear tooth under load at time t, δ mn (t) is the tooth profile deviation of each sheet gear at time t.

6. The method for calculating the meshing stiffness of a gear pair considering thermal-elastic coupling effect according to claim 5, characterized in that, The expression of the meshing stiffness is as follows: where j = 1, 2 denotes the driving and driven wheels; k tb is the face bending stiffness, k ts is the face shear stiffness, k ta is the face radial compression stiffness, k tf is the face radial matrix stiffness, k ab is the axial bending stiffness, k at is the axial torsional stiffness, k af is the axial matrix stiffness, k h is the Hertz contact stiffness, k is the mesh stiffness.

7. The method for calculating the meshing stiffness of a gear pair considering thermal-elastic coupling effect according to claim 1, characterized in that, S3 the friction coefficient μ of the tooth surface mn (t) calculated using load sharing theory and elastohydrodynamic lubrication theory; said engagement position relative sliding speed v mn (t) is calculated by the tooth surface slip rate formula.

8. The method for calculating the meshing stiffness of a gear pair considering thermal-elastic coupling effect according to claim 1, characterized in that, The expression of the comprehensive friction coefficient of the tooth surface is as follows: μ mn (t) = γ(t)μ a (t) + (1 - γ(t))μ b ; wherein μ a (t) is the oil film friction coefficient, calculated according to the elastohydrodynamic lubrication theory; μ b is the roughness peak friction coefficient, determined according to the material properties and roughness; t is time, μ mn (t) is the comprehensive friction coefficient at time t, and γ(t) is the oil film bearing ratio at time t.

9. The method for calculating the meshing stiffness of a gear pair considering thermal-elastic coupling effect according to claim 7, characterized in that, The calculation method of the oil film bearing ratio comprises: The initial oil film carrying ratio γ(t=0) = 1 is defined; iterative solution is carried out in combination with the thermal elastohydrodynamic lubrication analysis process until the rough peak contact pressure p c1 and p c2 The difference is less than the allowable value, that is, the actual oil film carrying ratio is obtained; where E0is the equivalent elastic modulus of the contact, R is the integrated radius of curvature of the contact point, L is the length of the line of tooth contact, β s is the roughness peak radius, n d is the roughness peak density, σ s , d d is the roughness height correlation number, where d d = 1.15σ s , h c is the central oil film thickness, W is the dimensionless load, F 5 / 2 is a function of the corresponding film thickness variation.

10. The method for calculating the meshing stiffness of a gear pair considering thermal-elastic coupling effect according to claim 5, characterized in that, In S4, the expression of the tooth surface thermal deformation is as follows: where x q is the x-direction thermal deformation of the gear tooth, y q is the y-direction thermal deformation of the gear tooth, z q is the z-direction thermal deformation of the gear tooth, T q is the gear node temperature, K S is the gear structure stiffness matrix, F TS is the equivalent node load matrix, [K S ] -1 is the inverse matrix of K S .

Images