Finite Element Mesh Generation Method for Helical Gears Based on MATLAB and ABAQUS

Abstract

The invention discloses a finite element mesh generation method for helical gears based on MATLAB and ABAQUS, comprising the following steps: 1) setting gear geometric parameters and helical gear finite element mesh size control parameters; 2) calculating the node coordinates on the involute of the end face of the helical gear based on the spatial position of the contact line; 3) simulating the gear hobbing process to obtain the envelope of the tooth root transition curve on the end face, and determining the node coordinates on the transition curve; 4) determining the node coordinates of different layer depths on the end face; 5) modifying the node coordinates of the end face surface according to the shaping method and shaping amount to obtain the modified node coordinates on the end face; 6) based on the modified node coordinates on the end face, obtaining the node coordinates on all cross sections, and then connecting these nodes into units: 7) obtaining the finite element mesh model of the helical gear. The invention can be used to study the influence of different tool parameters on tooth root stress, and is a powerful analysis tool for improving the bending strength of tooth root.

Description

Technical Field

[0001] The invention relates to the field of finite elements, and in particular to a finite element mesh generation method for helical gears based on MATLAB and ABAQUS. Background Art

[0002] Gear transmission has the advantages of compact structure, large load-bearing capacity, high efficiency and smooth movement. It is widely used in aerospace, metallurgical machinery, robotics, wind power, nuclear power and other new energy industries. The reliable operation of gear transmission has a significant impact on the economic benefits and life safety of the country and the people. It is of great significance to accurately analyze the stress state and reliability of gears. The industry mainly designs and strength checks gears based on international standards, national standards and industry standards. In the national standard GB / T3480.2-2021 (ISO6336-2:2019) "Calculation of load-bearing capacity of spur gears and helical gears Part 2: Calculation of tooth contact strength (pitting)", the contact pressure of the gear is used as the key parameter for tooth contact fatigue failure.

[0003] Meshing is a key step in finite element modeling, and the accuracy of finite element results depends on the quality of meshing. In finite element software, two methods can be used to generate mesh models. The first method is from geometry to mesh, which first generates the geometric model of the part, then partitions the part, lays seeds on the edge of the part, and generates a mesh. The geometric model can be established in finite element software or other 3D software. The second method is from node to mesh, which first generates nodes based on node information, and then connects the nodes into units to generate a mesh model. The first method relies on the function of geometric modeling software and is more suitable for simple parts. The self-geometry to mesh method can be used to generate standard spur gear and helical gear geometric models, and modified gear geometric models can also be generated, but these geometric models are often inaccurate and have poor mesh quality. In addition, if you want to perform special meshing, the first method is often difficult to achieve by combining coarse and fine meshes.

[0004] Compared with spur gears, the finite element model requirements of helical gears are higher. The line contact pressure of the helical gear finite element mesh using conventional partitioning methods is usually displayed as non-smooth multi-elliptical spots. Although the contact pressure can be made smoother by continuously refining the mesh size, this phenomenon of multi-elliptical spots cannot be fundamentally solved, and it increases storage and time costs. The parallelism of the edges of the nodes or units to the contact line is the basis for the accurate contact pressure distribution of helical gears. Conventional methods such as finite element-based meshing cannot achieve such a mesh distribution. It is necessary to develop a relatively general and accurate unconventional mesh generation method for helical gears. Summary of the invention

[0005] The purpose of the present invention is to provide a finite element mesh generation method for helical gears based on MATLAB and ABAQUS, comprising the following steps:

[0006] 1) Set the gear geometry parameters and the helical gear finite element mesh size control parameters;

[0007] 2) Based on the spatial position of the contact line, calculate the node coordinates on the involute of the helical gear end face;

[0008] 3) Simulate the gear hobbing process, obtain the envelope of the tooth root transition curve on the end face, and determine the node coordinates on the transition curve;

[0009] 4) Based on the node coordinates on the involute of the helical gear end face and the node coordinates on the transition curve, determine the node coordinates at different layer depths on the end face;

[0010] 5) modifying the node coordinates of the end surface layer according to the shaping method and the shaping amount to obtain the modified node coordinates on the end surface;

[0011] 6) Based on the modified node coordinates on the end face, use the lead spiral to scan the nodes on the end face to obtain the node coordinates on all cross sections, and then connect these nodes into units:

[0012] 7) According to the fixed format of the inp file, write the node number, coordinates, unit number and corresponding node number into the inp file;

[0013] 8) Import the inp file into ABAQUS to obtain the finite element mesh model of the helical gear.

[0014] Furthermore, the geometric parameters of the helical gear include but are not limited to module, number of teeth, pressure angle, helix angle, modification coefficient, rotation direction, and tool fillet radius.

[0015] The helical gear finite element mesh size control parameters include but are not limited to the axial size of the unit, the number of units in the near-surface area along the depth direction, the number of units in the far-from-surface area along the depth direction, and the number of units in the transition curve area along the tooth profile direction.

[0016] Further, the step of calculating the node coordinates on the involute of the helical gear end face includes:

[0017] 2-1) The intersection of the involute on the end face and the tooth top is taken as the first node, and then the coordinates of the second, third, ..., N nodes on the tooth top are determined according to the axial unit size; N is the number of tooth top nodes;

[0018] 2-2) Determine the coordinates of the jth node on the tooth profile through the intersection of the action surface of the jth node on the tooth top and the involute on the end face, until reaching the intersection of the involute on the end face and the transition curve; j = 2, 3..., M; M is the number of tooth profile nodes.

[0019] Further, the steps of determining the coordinates of the second, third, ..., Nth nodes on the tooth top include:

[0020] 2-1-1) Calculate the lead S of the helical gear, that is:

[0021]

[0022] Where: m t is the end face module, z is the number of teeth; β is the helix angle of the helical gear;

[0023] 2-1-2) According to the given unit axial dimension z e , calculate the angle α between the involutes of adjacent cross sections t ,Right now:

[0024] α t =z e / S (2)

[0025] 2-1-3) Calculate the pressure angle α corresponding to the i+1th node on the end face involute i+1 ,Right now:

[0026]

[0027] Where: α i+1 is the pressure angle corresponding to the i+1th node on the end face involute, r b is the base circle radius, r i is the radius of the i-th node on the end face involute; i = 1, 2, ..., N-1;

[0028] 2-1-4) Calculate the coordinates of the i+1th node (θ i+1 , r i+1 ),Right now:

[0029] The expansion angle corresponding to the i+1th node on the end face involute is:

[0030] θ i+1 =tan(α i+1 )-α i+1 (4)

[0031]

[0032] Among them, θ i+1 is the expansion angle corresponding to the i+1th node on the end face involute; r i+1 is the radius corresponding to the i+1th node.

[0033] Further, the steps of simulating the gear hobbing process, obtaining the envelope of the tooth root transition curve on the end face, and determining the node coordinates on the transition curve include:

[0034] 3-1) Discrete the involute extension line corresponding to the transition curve into n equal parts according to the radius, with a total of n+1 undetermined points;

[0035] 3-2) For the first node on the transition curve, the rack tool is moved horizontally. The entire rack tool translation process is discretized into m positions. The larger the m, the more continuous the process.

[0036] 3-3) According to the pure rolling relationship between the pitch line of the rack tool and the pitch circle of the gear, the gear is rotated accordingly, thereby causing the undetermined point and the generating line passing through the undetermined point to deflect, and recording the distance between the intersection of the generating line and the tool tooth profile and the undetermined point, a total of m distance values;

[0037] 3-4) Take the minimum value among these m distances, and the intersection point between the generation line corresponding to the minimum value and the tool tooth profile is the discrete point coordinate of the transition curve;

[0038] 3-5) Repeat the above process from 3-2) to 3-4) for the other n undetermined points, and discretize the transition curve into n+1 points;

[0039] 3-6) According to the number of transition curve nodes set, calculate the proportion of the arc length connecting each node to the first node, use the coordinates of n+1 discrete points to calculate the length of the transition curve, and calculate the arc length of each discrete point relative to the first discrete point and its proportion, compare the arc length proportions, and determine the node coordinates on the transition curve.

[0040] Further, in the process of horizontally moving the rack tool in step 3-2), the center of the gear is used as the coordinate origin, and the shape of the tool chamfer on the end face is an ellipse, wherein the major semi-axis a of the ellipse is in the rack translation direction, and the minor semi-axis b of the ellipse is perpendicular to the major semi-axis;

[0041] The equation of an ellipse is given by:

[0042]

[0043] Where: x tool and tool are the centers of the ellipse during tool movement; (x, y) are the coordinates of a point on the ellipse;

[0044] Ellipse major axis a, ellipse minor axis b, ellipse center y tool They are as follows:

[0045]

[0046] b=r tool (8)

[0047] y tool =rroot +r tool (9)

[0048] Where: r tool r is the tool corner radius; root is the gear tooth root circle radius.

[0049] Further, in step 3-3), the coordinates of the intersection of the line and the tool tooth profile (x p ,y p ) is as follows:

[0050]

[0051] y p =k i (x p -x i )+y i (11)

[0052] Among them, k i is the slope of the tangent line through the point; x i and i The coordinates of the unknown point of the i-th involute extension line;

[0053] The coefficient a1, coefficient b1, and coefficient c1 are as follows:

[0054]

[0055] Further, the step of determining the node coordinates of different layer depths on the end surface includes:

[0056] 4-1) Rotate the nodes in the involute area on the end face to obtain the node coordinates at different layer depths;

[0057] The node coordinates at different depths are rotated by angle α T and the rotation coordinates (x hi ,y hi ) characterization, namely:

[0058]

[0059] Where: x h0 and h0 are the coordinates of the nodes on the involute surface, x hi and hi is the coordinate of the depth node of the i-th layer in the involute region; h i is the depth of the i-th layer;

[0060] 4-2) According to the intersection of the first and last nodes of the transition curve and their normals, determine the center and radius of the transition curve. According to different depths, make concentric circles with different radii, translate the node coordinates of the array surface radially, and determine the node coordinates of different depths on the transition curve.

[0061] Further, in step 5), the modified node coordinates (x mi ,y mi ) is as follows:

[0062]

[0063] Among them, the angle Δα i As shown below:

[0064]

[0065] Among them, r b is the base circle radius;

[0066] Modification amount Δ i As shown below:

[0067]

[0068] Among them, m n is the normal modulus; L ca is the modification range coefficient; C aa is the maximum modification amount; r i is the radius corresponding to the i-th node;

[0069] The curvature radius pr corresponding to the starting point of the modification m0 As shown below:

[0070]

[0071] Where: r t is the radius of the tooth tip circle.

[0072]

[0073] Further, the step of connecting these nodes into units includes:

[0074] 6-1) The elements at the first and last end faces of the involute region are set as 6-node triangular prism wedge elements C3D6;

[0075] 6-2) The element in the middle of the involute region is set as an 8-node hexahedral reduced integration element C3D8R, and the face parallel to the tooth surface is a parallelogram;

[0076] 6-3) The unit in the tooth root transition curve area is set to an 8-node hexahedral reduced integration unit C3D8R, and the surface parallel to the tooth surface is a rectangle.

[0077] The technical effect of the present invention is undoubted. The beneficial effects of the present invention are:

[0078] 1. The present invention improves the mesh quality of the helical gear finite element mesh to provide a helical gear finite element model, which can obtain a more accurate, continuous and smooth contact pressure distribution;

[0079] 2. The present invention generates the finite element mesh of the helical gear by parameterization. By modifying the geometric parameters, the finite element mesh of the near tooth profile area of ​​the helical gear can be obtained within 30 seconds, saving a lot of pre-processing operations;

[0080] 3. The present invention includes a shaping method for linear tooth top trimming of helical gears, which can accurately obtain the finite element mesh of the helical gear after trimming, avoiding the process of three-dimensional geometric modeling and shaping in third-party software. The present invention can further expand other shaping methods such as drum shaping;

[0081] 4. The present invention obtains an accurate tooth root transition curve by simulating the process of gear hobbing, and can generate tooth root transition curves under different tool shapes, which can be further used to simulate and study the strength of different transition curves to improve the bending fatigue life of the tooth root;

[0082] 5. The present invention can set parameters for controlling the finite element mesh size of the helical gear in the program, which is convenient for quickly adjusting and generating finite element models of helical gears with different mesh sizes. Therefore, it can be reasonably set according to the configuration of the computer used and the required time.

[0083] 6. During the execution of the program in the present invention, some key intermediate results will be represented in the form of line graphs for rapid debugging of the program. BRIEF DESCRIPTION OF THE DRAWINGS

[0084] Figure 1 The overall flow chart of the method of the present invention is as follows;

[0085] Figure 2 is the node generation process on the end face, (2a) is the node on the tooth top, (2b) is the second node on the involute, and (2c) is multiple nodes on the involute;

[0086] Figure 3 Envelope formed for the gear profile and tool;

[0087] Figure 4 Generate the process for nodes on different cross sections.

[0088] Figure 5 It is the unit composition of the finite element mesh of the helical gear.

[0089] Figure 6 The contact pressure distribution of the gear finite element model generated based on the present invention. DETAILED DESCRIPTION

[0090] The present invention is further described below in conjunction with the embodiments, but it should not be understood that the above subject matter of the present invention is limited to the following embodiments. Without departing from the above technical ideas of the present invention, various substitutions and changes are made according to the common technical knowledge and customary means in the art, which should all be included in the protection scope of the present invention.

[0091] Embodiment 1:

[0092] See also Figures 1 to 6 , the finite element mesh generation method of helical gear based on MATLAB and ABAQUS includes the following steps:

[0093] 1) Set the gear geometry parameters and the helical gear finite element mesh size control parameters;

[0094] 2) Based on the spatial position of the contact line, calculate the node coordinates on the involute of the helical gear end face;

[0095] 3) Simulate the gear hobbing process, obtain the envelope of the tooth root transition curve on the end face, and determine the node coordinates on the transition curve;

[0096] 4) Based on the node coordinates on the involute of the helical gear end face and the node coordinates on the transition curve, determine the node coordinates at different layer depths on the end face;

[0097] 5) modifying the node coordinates of the end surface layer according to the shaping method and the shaping amount to obtain the modified node coordinates on the end surface;

[0098] 6) Based on the modified node coordinates on the end face, use the lead spiral to scan the nodes on the end face to obtain the node coordinates on all cross sections, and then connect these nodes into units:

[0099] 7) According to the fixed format of the inp file, write the node number, coordinates, unit number and corresponding node number into the inp file;

[0100] 8) Import the inp file into ABAQUS to obtain the finite element mesh model of the helical gear.

[0101] Embodiment 2:

[0102] The finite element mesh generation method for helical gears based on MATLAB and ABAQUS has the same technical content as Example 1. Furthermore, the geometric parameters of the helical gears include but are not limited to module, number of teeth, pressure angle, helix angle, displacement coefficient, rotation direction, and tool fillet radius.

[0103] The helical gear finite element mesh size control parameters include but are not limited to the axial size of the unit, the number of units in the near-surface area along the depth direction, the number of units in the far-from-surface area along the depth direction, and the number of units in the transition curve area along the tooth profile direction.

[0104] Embodiment 3:

[0105] The finite element mesh generation method for helical gears based on MATLAB and ABAQUS has the same technical content as any one of Embodiments 1-2. Further, the step of calculating the node coordinates on the involute of the end face of the helical gear includes:

[0106] 2-1) The intersection of the involute on the end face and the tooth top is taken as the first node, and then the coordinates of the second, third, ..., N nodes on the tooth top are determined according to the axial unit size; N is the number of tooth top nodes;

[0107] 2-2) Determine the coordinates of the jth node on the tooth profile through the intersection of the action surface of the jth node on the tooth top and the involute on the end face, until reaching the intersection of the involute on the end face and the transition curve; j = 2, 3..., M; M is the number of tooth profile nodes.

[0108] Embodiment 4:

[0109] The finite element mesh generation method for helical gears based on MATLAB and ABAQUS has the same technical content as any one of Embodiments 1 to 3. Further, the step of determining the coordinates of the second, third, ..., Nth nodes on the tooth top includes:

[0110] 2-1-1) Calculate the lead S of the helical gear, that is:

[0111]

[0112] Where: m t is the end face module, z is the number of teeth; β is the helix angle of the helical gear;

[0113] 2-1-2) According to the given unit axial dimension z e , calculate the angle α between the involutes of adjacent cross sections t ,Right now:

[0114] α t =z e / S (2)

[0115] 2-1-3) Calculate the pressure angle α corresponding to the i+1th node on the end face involute i+1 ,Right now:

[0116]

[0117] Where: α i+1is the pressure angle corresponding to the i+1th node on the end face involute, r b is the base circle radius, r i is the radius of the i-th node on the end face involute; i = 1, 2, ..., N-1;

[0118] 2-1-4) Calculate the coordinates of the i+1th node (θ i+1 , r i+1 ),Right now:

[0119] The expansion angle corresponding to the i+1th node on the end face involute is:

[0120] θ i+1 =tan(α i+1 )-α i+1 (4)

[0121]

[0122] Among them, θ i+1 is the expansion angle corresponding to the i+1th node on the end face involute; r i+1 is the radius corresponding to the i+1th node.

[0123] Embodiment 5:

[0124] The finite element mesh generation method for helical gears based on MATLAB and ABAQUS has the same technical content as any one of embodiments 1-4. Further, the steps of simulating the gear hobbing process, obtaining the envelope of the tooth root transition curve on the end face, and determining the node coordinates on the transition curve include:

[0125] 3-1) Discrete the involute extension line corresponding to the transition curve into n equal parts according to the radius, with a total of n+1 undetermined points;

[0126] 3-2) For the first node on the transition curve, the rack tool is moved horizontally. The entire rack tool translation process is discretized into m positions. The larger the m, the more continuous the process.

[0127] 3-3) According to the pure rolling relationship between the pitch line of the rack tool and the pitch circle of the gear, the gear is rotated accordingly, thereby causing the undetermined point and the generating line passing through the undetermined point to deflect, and recording the distance between the intersection of the generating line and the tool tooth profile and the undetermined point, a total of m distance values;

[0128] 3-4) Take the minimum value among these m distances, and the intersection point between the generation line corresponding to the minimum value and the tool tooth profile is the discrete point coordinate of the transition curve;

[0129] 3-5) Repeat the above process from 3-2) to 3-4) for the other n undetermined points, and discretize the transition curve into n+1 points;

[0130] 3-6) According to the number of transition curve nodes set, calculate the proportion of the arc length connecting each node to the first node, use the coordinates of n+1 discrete points to calculate the length of the transition curve, and calculate the arc length of each discrete point relative to the first discrete point and its proportion, compare the arc length proportions, and determine the node coordinates on the transition curve.

[0131] Embodiment 6:

[0132] A finite element mesh generation method for helical gears based on MATLAB and ABAQUS, the technical content of which is the same as any one of embodiments 1-5, further, in step 3-2), during the horizontal movement of the rack tool, the center of the gear is used as the coordinate origin, and the shape of the tool chamfer on the end face is an ellipse, wherein the major semi-axis a of the ellipse is in the rack translation direction, and the minor semi-axis b of the ellipse is perpendicular to the major semi-axis;

[0133] The equation of an ellipse is given by:

[0134]

[0135] Where: x tool and tool are the centers of the ellipse during tool movement; (x, y) are the coordinates of a point on the ellipse;

[0136] Ellipse major axis a, ellipse minor axis b, ellipse center y tool They are as follows:

[0137]

[0138] b=r tool (8)

[0139] y tool =r root +r tool (9)

[0140] Where: r tool r is the tool corner radius; root is the gear tooth root circle radius.

[0141] Embodiment 7:

[0142] The finite element mesh generation method for helical gears based on MATLAB and ABAQUS has the same technical content as any one of embodiments 1-6. Further, in step 3-3), the coordinates (x p ,y p ) is as follows:

[0143]

[0144] Among them, k iis the slope of the tangent line through the point; x i and i The coordinates of the unknown point of the i-th involute extension line;

[0145] The coefficient a1, coefficient b1, and coefficient c1 are as follows:

[0146]

[0147] Where: α r is the pressure angle of the end pitch circle.

[0148] Embodiment 8:

[0149] The finite element mesh generation method for helical gears based on MATLAB and ABAQUS has the same technical content as any one of Embodiments 1 to 7. Further, the step of determining the node coordinates of different layer depths on the end face includes:

[0150] 4-1) Rotate the nodes in the involute area on the end face to obtain the node coordinates at different layer depths;

[0151] The node coordinates at different depths are rotated by angle α T and the rotation coordinates (x hi ,y hi ) characterization, namely:

[0152]

[0153] Where: x h0 and h0 are the coordinates of the nodes on the involute surface, x hi and hi is the coordinate of the depth node of the i-th layer in the involute region; h i is the depth of the i-th layer;

[0154] 4-2) According to the intersection of the first and last nodes of the transition curve and their normals, determine the center and radius of the transition curve. According to different depths, make concentric circles with different radii, translate the node coordinates of the array surface radially, and determine the node coordinates of different depths on the transition curve.

[0155] Embodiment 9:

[0156] The finite element mesh generation method for helical gears based on MATLAB and ABAQUS has the same technical content as any one of embodiments 1-8. Further, in step 5), the modified node coordinates (x mi ,y mi ) is as follows:

[0157]

[0158] Among them, the angle Δα iAs shown below:

[0159]

[0160] Among them, r b is the base circle radius;

[0161] Modification amount Δ i As shown below:

[0162]

[0163] Among them, m n is the normal modulus; L ca is the modification range coefficient; C aa is the maximum modification amount; r i is the radius corresponding to the i-th node;

[0164] The curvature radius pr corresponding to the starting point of the modification m0 As shown below:

[0165]

[0166] Where: r t is the radius of the tooth tip circle.

[0167]

[0168] Embodiment 10:

[0169] The finite element mesh generation method for helical gears based on MATLAB and ABAQUS has the same technical content as any one of Embodiments 1 to 9. Further, the step of connecting these nodes into units includes:

[0170] 6-1) The elements at the first and last end faces of the involute region are set as 6-node triangular prism wedge elements C3D6;

[0171] 6-2) The element in the middle of the involute region is set as an 8-node hexahedral reduced integration element C3D8R, and the face parallel to the tooth surface is a parallelogram;

[0172] 6-3) The unit in the tooth root transition curve area is set to an 8-node hexahedral reduced integration unit C3D8R, and the surface parallel to the tooth surface is a rectangle.

[0173] Embodiment 11:

[0174] The finite element mesh generation method for helical gears based on MATLAB and ABAQUS is used to generate unconventional helical gear meshes to obtain more accurate contact pressure distribution, including the following steps:

[0175] 1) Set gear geometric parameters: Set the geometric parameters of the helical gear in the MATLAB program, including module, number of teeth, pressure angle, helix angle, displacement coefficient, rotation direction, tool fillet radius and other parameters;

[0176] 2) Setting the finite element mesh size control parameters of the helical gear: Setting the finite element mesh size control parameters in the MATLAB program, including the axial size of the unit, the number of units in the near-surface area along the depth direction, the number of units in the far-from-surface area along the depth direction, the number of units in the transition curve area along the tooth profile direction, and other parameters;

[0177] 3) Calculate the node coordinates on the end face involute based on the spatial position of the contact line:

[0178] 3-1) The intersection of the involute on the end face and the tooth top is taken as the first node;

[0179] 3-2) Determine the coordinates of the second, third, fourth, ... nodes on the tooth top according to the axial unit size;

[0180] 3-3) The coordinates of the second node on the tooth profile are determined by the intersection of the action surface of the second node on the tooth top and the involute on the end face, and the coordinates of the third node on the tooth profile are determined by the intersection of the action surface of the third node on the tooth top and the involute on the end face, and so on, until the intersection of the involute on the end face and the transition curve;

[0181] 4) Simulate the gear hobbing process, obtain the envelope of the tooth root transition curve on the end face, and determine the node coordinates on the transition curve:

[0182] 4-1) Discretize the involute extension line corresponding to the transition curve into n equal parts according to the radius, and set n to 200. This will allow the shape of the transition curve to be described more accurately, with a total of n+1 points to be determined;

[0183] 4-2) For the first node on the transition curve, the rack tool is moved horizontally. The entire rack tool translation process is discretized into m positions. The larger the m, the more continuous the process.

[0184] 4-3) According to the pure rolling relationship between the pitch line of the rack tool and the pitch circle of the gear, the gear is rotated accordingly, and the undetermined point and the generating line passing through the undetermined point are deflected, and the distance between the intersection of the generating line and the tool tooth profile and the undetermined point is recorded, with a total of m distance values;

[0185] 4-4) Take the minimum value among these m distances, and the intersection point between the corresponding generating line and the tool tooth profile is the coordinate of the discrete point of one of the transition curves;

[0186] 4-5) Repeat the above process from 4-2) to 4-4) for another n undetermined points, and the transition curve is discretized into n+1 points;

[0187] 4-6) According to the number of transition curve nodes set, calculate the proportion of the arc length connecting each node to the first node, use the coordinates of n+1 discrete points to calculate the length of the transition curve, and calculate the arc length of each discrete point relative to the first discrete point and its proportion, compare the arc length proportions, and determine the node coordinates on the transition curve.

[0188] 5) Determine the node coordinates of different depths on the end surface:

[0189] 5-1) Rotate the nodes in the involute area on the end face to obtain the node coordinates at different layer depths;

[0190] 5-2) According to the intersection of the first and last nodes of the transition curve and their normal lines, roughly determine the center and radius of the transition curve. According to different depths, make concentric circles with different radii, translate the node coordinates of the array surface radially, and determine the node coordinates of different depths on the transition curve;

[0191] 6) Modify the node coordinates according to the shaping method and shaping amount;

[0192] 7) Obtain the node coordinates on all cross sections by spirally scanning the nodes on the end surface according to the lead;

[0193] 8) Connect nodes into units:

[0194] 8-1) The elements at the first and last end faces of the involute region are set as 6-node triangular prism wedge elements (C3D6);

[0195] 8-2) The element in the middle of the involute region is set as an 8-node hexahedral reduced integration element (C3D8R), and the face parallel to the tooth surface is roughly a parallelogram;

[0196] 8-3) The element in the tooth root transition curve region is set as an 8-node hexahedral reduced integration element (C3D8R), and the surface parallel to the tooth surface is roughly rectangular;

[0197] 9) According to the fixed format of the inp file, write the node number, coordinates, unit number and corresponding node number into the inp file;

[0198] 10) Import the inp file into ABAQUS to obtain the finite element mesh model of the helical gear.

[0199] Step 3-2) When determining the coordinates of the 2nd, 3rd, 4th, ... nodes on the tooth top according to the axial unit size, the coordinates are determined by determining their radius, which is based on the helix angle β and lead S of the helical gear:

[0200]

[0201] Where: m tis the end module, and z is the number of teeth.

[0202] According to the given unit axial dimension z e , the angle between the involutes of adjacent cross sections is α t :

[0203] α t =z e / S (2)

[0204]

[0205] Where: α i+1 is the pressure angle corresponding to the i+1th node on the end face involute, r b is the base circle radius, r i is the radius of the i-th node on the end face involute. The expansion angle corresponding to the i+1-th node on the end face involute is:

[0206] θ i+1 =tan(α i+1 )-α i+1 (4)

[0207]

[0208] Step 4-2) During the horizontal movement of the rack tool, the center of the gear is used as the coordinate origin, and the chamfer of the tool is an ellipse on the end face, where the major semi-axis a is in the rack translation direction and the minor semi-axis b is perpendicular to the major semi-axis.

[0209]

[0210] b=r tool (7)

[0211] Where: r tool is the tool fillet radius. The corresponding ellipse equation is:

[0212]

[0213] Where: x tool and tool are the ellipse centers during tool movement, where x tool Controlled by the mobile process, y tool For a fixed value:

[0214] y tool =r root +r tool (9)

[0215] Where: r root is the gear tooth root circle radius.

[0216] In step 4-3), the distance between the intersection of the generating line and the tool tooth profile and the undetermined point is calculated. The intersection of the generating line and the ellipse conforms to a two-variable linear equation. The smaller of the two solutions is the correct result. The coefficients in the two-variable linear equation are as follows:

[0217]

[0218] Where: α r is the pressure angle of the end pitch circle.

[0219]

[0220] Where: x i and i The coordinates of the unknown point on the i-th involute extension line, k i is the slope of the tangent line through that point.

[0221]

[0222] The coordinates of the intersection of the tool tooth profile curve and the generating line passing through the point to be determined are as follows:

[0223]

[0224] y p =k i (x p -x i )+y i (14)

[0225] Step 5-1) Rotate the nodes in the involute area of ​​the end face to obtain the node coordinates of different layer depths, and calculate the rotation angle α T and rotated coordinates.

[0226]

[0227] Where: x h0 and h0 are the coordinates of the nodes on the involute surface, x hi and hi is the coordinate of the i-th depth node in the involute region.

[0228] Step 6) Modify the node coordinates according to the shaping method and shaping amount. The shaping method included in this method is linear trimming of the tooth top. The curvature radius corresponding to the starting point of the shaping is:

[0229]

[0230] Where: r t is the radius of the tooth tip circle, L ca is the modification range coefficient, m nis the normal modulus. The modification amount of different nodes is:

[0231]

[0232] The corresponding angles before and after the node modification are:

[0233]

[0234] The node coordinates after node modification are:

[0235]

[0236] Step 9) According to the fixed format of the inp file, write the node number, coordinates, unit number and corresponding node number into the inp file. The format of the inp file is as follows:

[0237] *Heading

[0238] **Job name:mesh_origin_2Model name:m0_global_mesh_or igin

[0239] **Generated by:Abaqus / CAE 2018

[0240] *Preprint,echo=NO,model=NO,history=NO,contact=NO

[0241] **

[0242] **PARTS

[0243] **

[0244] *Part,name=sun_surface_mesh

[0245] *Node

[0246] 101001,226.865835,12.269009,180.000000

[0247] 101002,226.244449,11.847940,180.000000

[0248] 101003,225.627589,11.435669,180.000000

[0249] 101004,225.015345,11.032117,180.000000

[0250] …

[0251] Node number, node x coordinate, node y coordinate, node z coordinate,

[0252] …

[0253] *Element,type=C3D6

[0254] 101001,301001,101002,101001,302001,102002,102001

[0255] 101002,301002,101003,101002,302002,102003,102002

[0256] 101003,301003,101004,101003,302003,102004,102003

[0257] …

[0258] Unit number, Node 1 number, Node 2 number, Node 3 number, Node 4 number, Node 5 number, Node 6 number

[0259] …

[0260] *Element,type=C3D8R

[0261] 101014,101014,301014,201015,101015,102014,302014,202015,102015

[0262] 102014,102014,302014,202015,102015,103014,303014,203015,103015

[0263] 103014,103014,303014,203015,103015,104014,304014,204015,104015

[0264] …

[0265] Unit number, Node 1 number, Node 2 number, Node 3 number, Node 4 number, Node 5 number, Node 6 number, Node 7 number, Node 8 number

[0266] …

[0267] *End Part

[0268] The name of the Job in the second row, the name of the Model, and the name of the Part in the eighth row can be modified. The node area is defined by the node number, followed by the x, y, and z coordinates of the node. This method sets up two types of units. The 6-node C3D6 area first defines the unit number, followed by the numbers of the 6 nodes that make up the unit. The definition of the unit in the 8-node C3D8R area is similar.

[0269] Step 10) Import the inp file into ABAQUS to obtain the finite element mesh model of the helical gear. The specific process is to open ABAQUS after generating the inp file, click File, Import, Part in the menu bar, and select the input inp file.

[0270] Embodiment 12:

[0271] This embodiment discloses a finite element mesh generation method for helical gears based on MATLAB and ABAQUS. This embodiment is a finite element mesh generation process for a helical gear. Figure 1 The process steps shown are carried out, which specifically include the following steps:

[0272] 1) Set the gear geometry parameters: Set the helical gear geometry parameters in the MATLAB program, as shown in Table 1:

[0273] Table 1 Geometric parameters of helical gears

[0274]

[0275] 2) Set the finite element mesh size control parameters of the helical gear: Set the finite element mesh size control parameters in the MATLAB program, as shown in Table 2:

[0276] Table 1 Grid size parameters of helical gears

[0277]

[0278] 3) Calculate the node coordinates on the end face involute based on the spatial position of the contact line:

[0279] 3-1) The intersection of the involute on the end face and the tooth top is taken as the first node;

[0280] 3-2) Determine the coordinates of the second, third, fourth, ... nodes on the tooth top according to the axial unit size;

[0281] Its coordinates are determined by determining its radius, based on the helix angle β and lead S of the helical gear:

[0282]

[0283] Where: m tis the end module, and z is the number of teeth.

[0284] According to the given unit axial dimension z e , the angle between the involutes of adjacent cross sections is α t :

[0285] α t =z e / S (2)

[0286]

[0287] Where: α i+1 is the pressure angle corresponding to the i+1th node on the end face involute, r b is the base circle radius, r i is the radius of the i-th node on the end face involute. The expansion angle corresponding to the i+1-th node on the end face involute is:

[0288] θ i+1 =tan(α i+1 )-α i+1 (4)

[0289]

[0290] 3-3) The coordinates of the second node on the tooth profile are determined by the intersection of the action surface of the second node on the tooth top and the involute on the end face. The coordinates of the third node on the tooth profile are determined by the intersection of the action surface of the third node on the tooth top and the involute on the end face. This process is as follows: Figure 2 As shown;

[0291] 4) Simulate the gear hobbing process to obtain the envelope of the tooth root transition curve on the end face, such as Figure 3 As shown, first determine the node coordinates on the transition curve:

[0292] 4-1) Discrete the involute extension line corresponding to the transition curve into 200 equal parts according to the radius, with a total of 201 undetermined points;

[0293] 4-2) For the first node on the transition curve, the rack tool is moved horizontally. The entire rack tool translation process is discretized into 10,000 positions. The involute is extended to intersect with the tooth root. This point is the end point of the tool processing. The starting point of the tool processing is the intersection of the gear tooth top circle and the involute.

[0294] The center of the gear is taken as the coordinate origin, and the tool chamfer is an ellipse on the end face, in which the major semi-axis a is in the rack translation direction and the minor semi-axis b is perpendicular to the major semi-axis.

[0295]

[0296] Where: r tool is the tool fillet radius. The corresponding ellipse equation is:

[0297]

[0298] Where: x tool and tool are the ellipse centers during tool movement, where x tool Controlled by the mobile process, y tool For a fixed value:

[0299] y tool =r root +r tool (9)

[0300] Where: r root is the gear tooth root circle radius.

[0301] 4-3) According to the pure rolling relationship between the pitch line of the rack tool and the pitch circle of the gear, the gear is rotated accordingly, and the pending point and the generating line passing through the pending point are deflected. The distance between the intersection of the generating line and the tool tooth profile and the pending point is recorded, with a total of 10,000 distance values. The intersection of the generating line and the ellipse conforms to a two-variable linear equation. The smaller of the two solutions is the correct result. The coefficients in the two-variable linear equation are as follows:

[0302]

[0303] Where: α r is the pressure angle of the end pitch circle.

[0304]

[0305] Where: x i and i The coordinates of the unknown point on the i-th involute extension line, k i is the slope of the tangent line through that point.

[0306]

[0307] The coordinates of the intersection of the tool tooth profile curve and the generating line passing through the point to be determined are as follows:

[0308]

[0309] y p =k i (p x -x i )+y i (14)

[0310] 4-4) Take the minimum value among these 10,000 distances, and the intersection point between the corresponding generating line and the tool tooth profile is the coordinate of a discrete point of one of the transition curves;

[0311] 4-5) Repeat the process from 4-2) to 4-2) for another 200 undetermined points, and the transition curve is discretized into 201 points;

[0312] 4-6) According to the number of transition curve nodes set, the arc length ratio of each node to the first node is calculated, and the length of the transition curve is calculated using the coordinates of 201 discrete points. The arc length of each discrete point relative to the first discrete point and its ratio are calculated. The arc length ratios are compared to determine the node coordinates on the transition curve.

[0313] 5) Determine the node coordinates of different depths on the end surface:

[0314] 5-1) Rotate the nodes in the involute area on the end face to obtain the node coordinates of different layer depths, and calculate the rotation angle α T and rotated coordinates.

[0315]

[0316]

[0317] Where: x h0 and h0 are the coordinates of the nodes on the involute surface, x hi and hi is the coordinate of the i-th depth node in the involute region.

[0318] 5-2) According to the intersection of the first and last nodes of the transition curve and their normal lines, roughly determine the center and radius of the transition curve. According to different depths, make concentric circles with different radii, translate the node coordinates of the array surface radially, and determine the node coordinates of different depths on the transition curve;

[0319] 6) Modify the node coordinates according to the shaping method and shaping amount. The shaping method included in the present invention is linear trimming of the tooth top, and the curvature radius corresponding to the starting point of the shaping is:

[0320]

[0321] Where: r t is the radius of the tooth tip circle, L ca is the modification range coefficient, m n is the normal modulus. The modification amount of different nodes is:

[0322]

[0323] The corresponding angles before and after the node modification are:

[0324]

[0325] The node coordinates after node modification are:

[0326]

[0327] 7) According to the lead spiral scanning nodes on the end surface, the node coordinates on all cross sections are obtained, such as Figure 4 As shown;

[0328] 8) Connect nodes into units:

[0329] 8-1) The elements at the first and last end faces of the involute region are set as 6-node triangular prism wedge elements (C3D6);

[0330] 8-2) The element in the middle of the involute region is set as an 8-node hexahedral reduced integration element (C3D8R), and the face parallel to the tooth surface is roughly a parallelogram;

[0331] 8-3) The unit in the tooth root transition curve area is set as an 8-node hexahedral reduced integration unit (C3D8R), and the surface parallel to the tooth surface is roughly rectangular, such as Figure 5 As shown;

[0332] 9) According to the fixed format of the inp file, write the node number, coordinates, unit number and corresponding node number into the inp file. Some inp files are as follows:

[0333] *Heading

[0334] **Job name:mesh_origin_2Model name:m0_global_mesh_or igin

[0335] **Generated by:Abaqus / CAE 2018

[0336] *Preprint,echo=NO,model=NO,history=NO,contact=NO

[0337] **

[0338] **PARTS

[0339] **

[0340] *Part,name=sun_surface_mesh

[0341] *Node

[0342] 101001,226.865835,12.269009,180.000000

[0343] 101002,226.244449,11.847940,180.000000

[0344] 101003,225.627589,11.435669,180.000000

[0345] 101004,225.015345,11.032117,180.000000

[0346] …

[0347] *Element,type=C3D6

[0348] 101001,301001,101002,101001,302001,102002,102001

[0349] 101002,301002,101003,101002,302002,102003,102002

[0350] 101003,301003,101004,101003,302003,102004,102003

[0351] …

[0352] *Element,type=C3D8R

[0353] 101014,101014,301014,201015,101015,102014,302014,202015,102015

[0354] 102014,102014,302014,202015,102015,103014,303014,203015,103015

[0355] 103014,103014,303014,203015,103015,104014,304014,204015,104015

[0356] …

[0357] *End Part

[0358] 10) Import the inp file into ABAQUS to obtain the finite element mesh model of the helical gear. The specific process is to open ABAQUS after generating the inp file, click File, Import, Part in the menu bar, and select the input inp file.

[0359] The contact pressure distribution of the gears at a certain meshing moment is as follows: Figure 6 As shown, the contact pressure is continuous, smooth and has high precision.

[0360] Embodiment 13:

[0361] The main steps of this embodiment are the same as those of Embodiment 12. Furthermore, the geometric parameters such as the module, number of teeth, pressure angle, displacement coefficient, tooth width and hand direction in step 1) are modified.

[0362] Embodiment 14:

[0363] The main steps of this embodiment are the same as those of Embodiment 12. Furthermore, the mesh control parameters in step 2) such as the axial size of the unit, the number of units in the depth direction, and the number of units in the tooth root transition curve are modified.

[0364] Embodiment 15:

[0365] The main steps of this embodiment are the same as those of Embodiment 12. Furthermore, the shaping parameters such as the maximum shaping amount and shaping range coefficient in step 6) are modified.

Claims

1. Finite element mesh generation method for helical gears based on MATLAB and ABAQUS, characterized by , including the following steps: 1) Set the gear geometry parameters and the helical gear finite element mesh size control parameters; 2) Based on the spatial position of the contact line, calculate the node coordinates on the involute of the helical gear end face; 3) Simulate the gear hobbing process, obtain the envelope of the tooth root transition curve on the end face, and determine the node coordinates on the transition curve; 4) Based on the node coordinates on the involute of the helical gear end face and the node coordinates on the transition curve, determine the node coordinates at different layer depths on the end face; 5) Modifying the node coordinates on the end surface according to the shaping method and shaping amount to obtain the modified node coordinates on the end surface; 6) Based on the modified node coordinates on the end face, use the lead spiral to scan the nodes on the end face to obtain the node coordinates on all cross sections, and then connect these nodes into units: 7) According to the fixed format of the inp file, write the node number, coordinates, unit number and corresponding node number into the inp file; 8) Import the inp file into ABAQUS to obtain the finite element mesh model of the helical gear; The steps of simulating the gear hobbing process, obtaining the envelope of the tooth root transition curve on the end face, and determining the node coordinates on the transition curve include: 3-1) Discrete the involute extension line corresponding to the transition curve into n equal parts according to the radius, with a total of n+1 undetermined points; 3-2) For the first node on the transition curve, the rack tool is moved horizontally, and the entire rack tool translation process is discretized into m positions; 3-3) According to the pure rolling relationship between the pitch line of the rack tool and the pitch circle of the gear, the gear is rotated accordingly, thereby causing the undetermined point and the generating line passing through the undetermined point to deflect, and recording the distance between the intersection of the generating line and the tool tooth profile and the undetermined point, a total of m distance values; 3-4) Take the minimum value among these m distances, and the intersection point between the generation line corresponding to the minimum value and the tool tooth profile is the discrete point coordinate of the transition curve; 3-5) Repeat the above process from 3-2) to 3-4) for the other n undetermined points, and discretize the transition curve into n+1 points; 3-6) According to the number of transition curve nodes set, calculate the proportion of the arc length connecting each node to the first node, use the coordinates of n+1 discrete points to calculate the length of the transition curve, and calculate the arc length of each discrete point relative to the first discrete point and its proportion, compare the arc length proportions, and determine the node coordinates on the transition curve.

2. The finite element mesh generation method for helical gears based on MATLAB and ABAQUS according to claim 1 is characterized in that: The geometric parameters of the helical gear include module, number of teeth, pressure angle, helix angle, modification coefficient, rotation direction, and tool fillet radius; The helical gear finite element mesh size control parameters include the axial size of the unit, the number of units in the near-surface area along the depth direction, the number of units in the far-from-surface area along the depth direction, and the number of units in the transition curve area along the tooth profile direction.

3. The finite element mesh generation method for helical gears based on MATLAB and ABAQUS according to claim 1 is characterized in that: The steps for calculating the node coordinates on the involute of the helical gear end face include: 2-1) The intersection of the involute on the end face and the tooth top is taken as the first node, and then the coordinates of the second, third, ..., N nodes on the tooth top are determined according to the axial unit size; N is the number of tooth top nodes; 2-2) Determine the coordinates of the jth node on the tooth profile through the intersection of the action surface of the jth node on the tooth top and the involute on the end face, until reaching the intersection of the involute on the end face and the transition curve; j = 2, 3..., M; M is the number of tooth profile nodes.

4. The finite element mesh generation method for helical gears based on MATLAB and ABAQUS according to claim 3 is characterized in that: The steps of determining the coordinates of the second, third, ..., Nth nodes on the tooth top include: 2-1-1) Calculate the lead S of the helical gear, that is: Where: m t is the end face module, z is the number of teeth; β is the helix angle of the helical gear; 2-1-2) According to the given unit axial dimension z e , calculate the angle α between the involutes of adjacent cross sections t ,Right now: α t =z e / S (2) 2-1-3) Calculate the pressure angle α corresponding to the i+1th node on the end face involute i+1 ,Right now: Where: α i+1 is the pressure angle corresponding to the i+1th node on the end face involute, r b is the base circle radius, r i is the radius of the i-th node on the end face involute; i = 1, 2, ..., N-1; 2-1-4) Calculate the coordinates of the i+1th node (θ i+1 , r i+1 ),Right now: The expansion angle corresponding to the i+1th node on the end face involute is: i i+1 =tan(α i+1 )-a i+1 (4) Among them, θ i+1 is the expansion angle corresponding to the i+1th node on the end face involute; r i+1 is the radius corresponding to the i+1th node.

5. The finite element mesh generation method for helical gears based on MATLAB and ABAQUS according to claim 1 is characterized in that: In the process of horizontally moving the rack tool in step 3-2), the center of the gear is used as the coordinate origin, and the shape of the tool chamfer on the end face is an ellipse, wherein the major semi-axis a of the ellipse is in the rack translation direction, and the minor semi-axis b of the ellipse is perpendicular to the major semi-axis; The equation of an ellipse is given by: Where: x tool and tool are the centers of the ellipse during tool movement; (x, y) are the coordinates of the points on the ellipse; Ellipse major axis a, ellipse minor axis b, ellipse center y tool They are as follows: b＝r tool (8) y tool ＝r root +r tool (9) Where: r tool r is the tool corner radius; root is the gear tooth root circle radius.

6. The finite element mesh generation method for helical gears based on MATLAB and ABAQUS according to claim 5 is characterized in that: In step 3-3), the coordinates of the intersection of the generating line and the tool tooth profile (x p ,y p ) is as follows: and p =k i (x p -x i )+y i (11) Among them, k i is the slope of the tangent line through the point; x i and i The coordinates of the unknown point of the i-th involute extension line; The coefficient a1, coefficient b1, and coefficient c1 are as follows:

7. The finite element mesh generation method for helical gears based on MATLAB and ABAQUS according to claim 1 is characterized in that: The steps for determining the node coordinates at different depths on the end surface include: 4-1) Rotate the nodes in the involute area on the end face to obtain the node coordinates at different layer depths; The node coordinates at different depths are rotated by angle α T and the rotation coordinates (x hi ,y hi ) characterization, namely: Where: x h0 and h0 are the coordinates of the nodes on the involute surface, x hi and hi is the coordinate of the depth node of the i-th layer in the involute region; h i is the depth of the i-th layer; 4-2) According to the intersection of the first and last nodes of the transition curve and their normals, determine the center and radius of the transition curve. According to different depths, make concentric circles with different radii, translate the node coordinates of the array surface radially, and determine the node coordinates of different depths on the transition curve.

8. The finite element mesh generation method for helical gears based on MATLAB and ABAQUS according to claim 1 is characterized in that: In step 5), the modified node coordinates (x mi ,y mi ) is as follows: Among them, the angle Δα i As shown below: Among them, r b is the base circle radius; Modification amount Δ i As shown below: Among them, m n is the normal modulus; L ca is the modification range coefficient; C aa is the maximum modification amount; r i is the radius corresponding to the i-th node; The curvature radius pr corresponding to the starting point of the shaping m0 As shown below: Where: r t is the radius of the tooth tip circle.

9. The finite element mesh generation method for helical gears based on MATLAB and ABAQUS according to claim 1, characterized in that: The steps to connect these nodes into units include: 6-1) The elements at the first and last end faces of the involute region are set as 6-node triangular prism wedge elements (C3D6); 6-2) The element in the middle of the involute region is set as an 8-node hexahedral reduced integration element (C3D8R), and the face parallel to the tooth surface is a parallelogram; 6-3) The element in the tooth root transition curve area is set as an 8-node hexahedral reduced integration element (C3D8R), and the surface parallel to the tooth surface is a rectangle.

Images