Abstract
In order to perform a tooth contact analysis of helical gears with satisfactory accuracy and computational time consuming, a parameterized approach to establish a high-precision three-dimension (3D) finite element model (FEM) of involute helical gears is proposed. The enveloping theory and dentiform normal method are applied to deduce the mathematical representations of the root transit curve as well as the tooth profile of the external gear in the transverse plane based on the manufacturing process. A bottom-up modelling method is applied to build the FEM of the helical gear directly without the intervention of CAD software or creating the geometry model in advance. Local refinement methodology of the hexahedral element has been developed to improve the mesh quality and accuracy. A computer program is developed to establish 3D helical gear FEM with contact region well refined with any parameters and mesh density automatically. The comparison of tooth contact analysis between the coarse-mesh model and local refinement model demonstrates that the present method can efficiently improve the simulation accuracy while greatly reduce the computing cost. Using the proposed model, the tooth load sharing ratio, static transmission error, meshing stiffness, root bending stress, and contact stress of the helical gear are obtained based on the quasistatic load tooth contact analysis. This methodology can also be used to create other types of involute gears, such as high contact ratio gear, involute helical gears with crossed axes, or spiral bevel gears.
1. Introduction
Gear-loaded tooth contact analysis (LTCA) is important for the design and analysis of gear performance. Many researchers have proposed several methods to predict the contact pattern and the gear mesh compliance: experimental method, such as photo elastic measurement and ink marking [1, 2], analytical methods [3–6], applying the elastic contact theory for LTCA, and finite element method [7–10]. From all the methods, FEM is the most effective calculation technique to predict the root bending stress, contact stress distribution, static transmission error, and mesh stiffness [11–15], neglecting the difficulties of preparing FEM beforehand. A two-dimensional (2D) spur gear FEM is built by Wang [16] to calculate the static transmission error, single and combined torsional mesh stiffness, tooth load-sharing ratio, maximum tooth root stress, and surface contact stress against various input torques over a complete mesh cycle. However, there is a large difference in the meshing stiffness calculation of the helical gear from the spur gear. Since the contact line of a helical gear is inclined, the load distribution on the contact line is nonuniform. Therefore, the LTCA of the helical gear must be modeled in three dimensions (3D), leading to more complicated calculations.
The modelling accuracy and element size have an important impact on the results of LTCA. In order to perform an LTCA with satisfactory accuracy, the mesh grid must be very fine that there would be considerable difficulties to solve the case. These problems make the contact analysis of helical gears a very challenging task. The unique features of the tooth surface contact make the traditional finite element method rather inefficient and inaccurate. For gear contact analysis, the mesh size in the contact area must be far smaller than the contact width, which requires a very fine mesh near the contact area. When performing contact analysis for the entire meshing cycle, well-refined mesh is indispensable over the whole tooth surfaces. The computational resource requirements and calculation time required to accommodate such fine meshed model are absolutely unacceptable. In addition, the traditional modelling method is far from the modelling accuracy required for tooth surface contact analysis. Moreover, traditional finite element methods have difficulties in generating controllable and optimal 3D meshes that is capable of modelling the stress gradients in critical areas while minimizing the number of total elements. In order to perform an LTCA with satisfactory accuracy while reducing the computational cost, a high-precision 3D FEM with local contact mesh refinement is necessary.
Extensive studies have focused on the LTCA of helical gears with FEM, and significant progresses have been made. Qin [17] and Gurumani and Shanmugam [10] developed a 3D helical gear FEM with the addition of computer-aided design (CAD) software and imported it into a finite element solver for static simulation. Contact load distribution and contact stress are also calculated to analyze the crowning effects. The CAD/FEA method of geometry generation causes the loss of geometry accuracy and the modelling parameters. It is acceptable for bending calculations, but for contact calculations, the tooth flanks are often not accurate enough. In order to get high-accuracy gear surface, Mao [18] generated the gear tooth profile node coordinates in C++ according to the gear meshing theory, imported the data into the finite element analysis software through Python script to generate the 2D gear geometric model, and then extruded the plane model to obtain a 3D helical gear geometry to investigate the gear surface fatigue wear phenomena. The method is only suitable for spur gears or occasions where the accuracy requirement is not very strict. Since the extrude method is unreliable in finite element software, edge contact may occur in 3D helical gear contact analysis when the helix angle or face width is large. Further development has been made by Hedlund [19] who utilized eight-node trilinear hexahedral elements that stretched automatically over a determined domain to generate tooth surface geometry. Attempts have been made by Brauer [20], Patil [21], and Jao Huang [22] to create the FE model through meshing the solid model, which was built from the bottom-up and divided into several regular quadrilateral regions beforehand. In these methods, the geometric model had been generated with commercial finite element software directly instead of delivering from other CAD packages. But curve fitting or interpolation is still required in mesh grid generation and may cause deviation of the nodes on the tooth surface. Recently, Barbieri [23] obtained the tooth profile by numerical simulation of the cutting process and provided an accurate description in NURBS curves to generate a 3D helical gear model. Refinement criteria are applied in 2D FEM in order to obtain accurate solutions in LTCA. However, the coarse-mesh model is still used for contact analysis of 3D FEM. In many other research publications [24–26], a coarse meshed finite element gear model is also used for gear contact simulation. In order to make the result of contact analysis more accurate, the mesh density of the entire tooth must be increased [8, 27, 28]. The main disadvantage of these methods is the great computing resource requirement and dramatic increase in computation time when the number of elements is very large.
Although some progresses have been made in modelling a relatively high-precision 3D helical gear finite model with different methods, most of the models are coarse mesh and none of these studies make further research on the local contact pattern due to the difficulties in hexahedral mesh refinement on the tooth contact surface. Therefore, further studies are still essential.
The main objective in this paper is to present a method to establish parameterized 3D FEM of involute spur/helical gears with high accuracy. The mathematical equations of tooth profiles of the external gear would be obtained based on the profile of the machine tool, the manufacturing process, and gear meshing theory. Discrete nodes of one tooth will be generated on the tooth profile and transition curve in the transverse section in ANSYS. Eight-node hexahedral element solid 185 is applied to build the element model directly, without creating a geometric model firstly or generating the tooth surface nodes with the interaction of other software, i.e., C++ or MATLAB. On the basis of the initial precision coarse-mesh model, an efficiently localized refinement of the elements in the contact region on the bilateral flanks will be introduced to ensure the mesh quality for contact analysis. An advanced nonlinear tooth contact analysis will be applied to investigate gear contact behaviors based on the helical gear FEM proposed in the article.
2. Tooth Profile Equations in the Transverse Plane
An accurate tooth profile on the transverse plane is necessary to establish a high-precision finite element helical gear model. When calculating the transmission error and mesh stiffness of the spur and helical gear, the precision of the root fillet curve has great influence on the tooth deflection. Different from the working profile, which can be easily established by using the standard involute equation, the root transition curve can only be generated by the fillet curve of the rack cutter. In this paper, the enveloping theory and dentiform normal method are applied to deduce the equations of the root transit curve as well as the involute tooth profile of the gear.
2.1. Transverse Tooth Profile of the External Gear
Figure 1(a) shows the tooth profile parameters in the transverse section of a rack cutter. It is formed by using three parts in one side, the oblique line (part III) generates the involute profile of the gear, the horizontal line (part I) generates the dedendum circle of the gear, and the arc of radius (part II) connects part I and part III. In the rigidly connected coordinate system of the rack cutter, the fillet curve is represented by the following equations:wherewhere is the transverse module, is the transverse pressure angle, and are the transverse addendum coefficient and tip clearance coefficient, and is the transverse modification coefficient.
(a)
(b)
Part III is represented by the following equation:
The superscript indicates that the values are measured in the transverse section. These equations are equally applicable to the spur gear.
The coordinate systems , , and shown in Figure 1(b) are used to derive the tooth profile equations. The coordinate systems and are connected to the rack cutter and the gear to perform translational and rotational motions with respect to the fixed coordinate system. The coordinate transformation from to is based on the following matrix equation:and the translational matrix of the rack gear is
The generated profile of the gear in the coordinate is based on the theorem that a point of the rack cutter generates the respective point of the gear at a position where the normal to the shape of the rack cutter passes through the instantaneous center of rotation [29]. Therefore, the displacement of the rack cutter and the point on the rack cutter have the relationship as
In the fillet curve part (arc II) of the rack cutter, the tangent angle of the arc and its parameter are related bywhere is the angle between the tangent line of the arc at point with the x-axis.
Substituting (1), (2), and (7) in (6) yields the following equation of the contact line represented by the arc parameter :
According to Figure 1, the rotation angle of the gear and the displacement of the rack cutter are related by , and then the meshing equation of the rack cutter and the gear is represented by
Then, the equation of the transit curve of the gear represented by the arc parameter and gear rotation angle is
As to the oblique line (part III),
Substituting (3) and (11) into (6) yields
Then, the equation of the involute of gear represented by using the parameter of the rack cutter and rotation angle is given by
3. Finite Element Modelling Methodology of 3D Helical Gear
The parameters of the gears are given in Table 1.
The flowchart shown in Figure 2 gives the parameterized finite element modelling method for the 3D external helical gear.
Quadrilateral and hexahedron meshes are the most ideal grid structures in finite element structural analysis. Hexahedral meshes have incomparable advantages over tetrahedral meshes in terms of computational accuracy, number of elements and nodes generated by using hexahedral mesh, memory requirement and computation time, and deformation characteristics. Therefore, the three-dimensional FEM of helical gears is the hexahedral mesh.
The tooth profile is established based on the equations deduced in Section 2 as shown in Figure 3. To facilitate parameterized control of the nodes and the mesh size and shape and use the high-quality 3D element, the tooth profile of the helical gear in the transverse section is divided into 8 regular quadrilateral area regions. Points and are the tangential points of the root transition curve and the dedendum circle of both sides, respectively, points ② and ⑤ are the tangential points of the root transition curve and the involute tooth profile, and points ① and ④ are the midpoint of the adjacent tooth root circle.
According to the requirement of the analysis, the number of keypoints (KPs) filled on the root fillet curve between KP① and KP② is supposed to be , on the involute contact profile between KP② and KP③ is , on the addendum circle between KP③ and KP⑨ is , and on the web sideline between KP① and KP⑫ is . Keypoint KP① is set at the center rotation of the gear.
According to the distribution of the KPs in Figure 3 and the number of elements to be created on the corresponding lines, the number of KP①–KP⑮ is presented as
Figure 4(a) shows the KPs filled on the tooth profile. To include the effect of the Hertzian deformation, the number of KPs created on the involute curve is based on the recommendation given by Coy [30]:where is the arc length of the involute curve.where and are the length and width of the element, respectively, and is the minimum Hertzian contact width in one mesh cycle:
(a)
(b)
(c)
(d)
(e)
(f)
A loop program is used to identify the number of KPs on the midline of the gear tooth corresponding to the KPs on the tooth profile and fill KPs between them as shown in Figure 4(b). Because of the high stress gradient of contact on the surface, the KPs are set to be from dense to sparse in the tooth thickness direction by controlling the proportional factor of the fill command (KFILL). The total number of KPs in the transverse plane is
Nodes are generated at every KPs on the end face. According to the direction of rotation and the helix angle, the nodes are copied-translated-rotated to generate the multilayer node model (Figure 4(c)). In order to guarantee the aspect ratio of the element on the flank, the layers of nodes in the face-width direction are determined aswhere is the face width of the gear.
An hexahedral element is defined for every 8 neighbouring nodes to generate the FEM of one tooth (Figure 4(d)). Then the complete 3D helical gear FEM can be established by copying the single-tooth FEM and rotate it around the axis as shown in Figure 4(e).
The mating gear is established in the modelling coordinate system and set up according to the center distance and the installation error\eccentricity error. The 3D helical gear pair FEM is shown in Figure 4(f).
4. Contact Region Refinement Technology
According to Hertz theory, an elliptical contact area is formed on the tooth surface when a pair of helical gears contacts each other. Although the recommended grid size has been given by Coy in [30], further studies are required to be conducted. On the one hand, his work did not include the grid size selection of three-dimensional contact problem; on the other hand, although the whole gear three-dimensional FEM build with the recommended grid size is very large, there are no more than two elements in the direction of contact width, and this mesh density is not fine enough to investigate the contact deformations of the tooth. In order to analyze the load distribution and contact pattern with considerable accuracy, the elements in the local contact area must be well refined.
4.1. Normalization of the Contact Region Elements
The coordinate system , whose y-axis is parallel to the line of action (LOA), is established at the center of the gear, as shown in Figure 5. Contacts will happen when the working side of the tooth is within the actual line of action . The nodes of each layer on the working side of the tooth flank nearest to the face of action (FOA) are selected and moved to the FOA by modifying the x coordinate value as . The numbering scheme of the nodes on the involute curve, namely, the node number increases with the radius of the node, makes it easy to select the ambilateral node in the same node layer, i.e., and . All the nodes are moved to the sides of the FOA separately to form a standard contact region that is parallel to the FOA as shown in Figure 6.
The elements of the contact region are classified into 3 types: type-I element group (EG)—only one element beside the FOA in one layer (the red element in Figure 6(d)); type-II EG—there are 2 elements beside the FOA (the green elements in Figure 6(d)); and type-III EG—there are 3 elements beside the FOA (the blue elements in Figure 6(d)). All of the element groups (EGs) may appear depending on the helix angle and the relative position of the flank and the FOA.
A bottom-up modelling method of the 3D helical gear FEM makes all the nodes parameterized. Intermediate data files are used to store the node numbers of the EGs in the numbering scheme shown in Figure 7. Once we get the minimum node number and give it to , the number of other nodes in the remaining corner of the EG can be represented as
(a)
(b)
(c)
Only two types of EGs will appear in spur gear transmission, type-I on the tooth tip or near the tooth root and type-II on the other parts except for the root and the tip. For the helical gear, type-III EG will definitely appear as transition connection of two type-II EGs. Three coarse elements are deleted, and two degenerate triangular prism elements are reconstructed along the diagonal line at the top and bottom of the contact zone, respectively, as shown in Figure 8. The rest of the middle part is rebuilt into a type-II EG. After doing this, all the EGs in the contact zone are reformed and only composed of type-II EGs.
(a)
(b)
4.2. Hexahedral Element Refinement Template
The ANSYS program allows to refine the mesh locally around specified nodes, elements, keypoints, lines, or areas for volume meshes composed of tetrahedra. Meshes composed of volume elements other than tetrahedra (for example, hexahedra, wedges, and pyramids) cannot be locally refined [31]. Therefore, it is necessary to establish reasonable transition element manually to refine the hexahedral mesh in the contact area. The retain element types method is applied to convert the elements into smaller elements of the same type.
A general parameterized hexahedral mesh refinement template with well-shaped and uniform-sized subdivided elements is established according to the refinement level as shown in Table 2. Refinement level 1 is the local node-refinement template and is not the focus of this article. For the refinement level not less than two, edge-refinement templates are developed at the edges coinciding with the contact line. The relationships between the number of nodes filled on the edge along the contact line and the number of subelements with the refinement level are
This refinement process is repeated along the contact line until all the original elements in the contact region are well refined, as is shown in Figure 9. It can be seen from the figure clearly that the elements in the contact region become fine and regular after refinement not only on the surface but also with a certain depth. That is crucial for contact analysis since the contact stress gradient near the surface would be very large.
4.3. Mapping the Surface Nodes to the Theoretical Position
In the refinement process, the nodes of the subelements are generated by linear interpolation on the edges of the coarse mesh. At the same time, the nodes on the tooth surface deviate from the theoretical position due to the movement of the nodes near the LOA to the contact strip. However, for tooth contact analysis, small geometric error can cause huge differences in results. Therefore, it is necessary to correct the coordinates of the nodes to the theoretical position of tooth surface.
Figure 10 illustrates the coordinate system for node mapping. A cylindrical coordinate system is set up at the end of the helical gear, and the x-axis is the symmetry axis of the tooth profile at the end face. The cylindrical coordinate system can be gotten by translating along the z-axis of distance and rotating by angle . Since -axis of the new system is still the symmetry axis of the tooth profile, the relationship between the rotation angle and the helix angle iswhere is the z coordinate value of the node on the tooth surface in , is the pitch radius, and is the number of teeth from the reference position.
Point in Figure 11(a) is a node generated by interpolation which is supposed to be on the involute tooth profile of the tooth surface. In order to mapping it to the theoretical position, geometric relationship of the interpolation point and theoretical point of the external gear is shown in Figure 11(b). The angle between line and the horizontal -axis can be obtained through the trigonometric function by extracting the polar radius and vertical coordinate value of the interpolation point . Arc mapping point is generated by the intersection of the arc with radius and the theoretical involute. Since point is on the theoretical tooth profile, the theoretical pressure angle, which has the relationship , and can also be obtained by applying gear geometry theory. The offset angle between the arc mapping point and point is , . The central angle of the node with a polar radius of is calculated based on the theoretical tooth thickness formula of the gear. The angle between the generating line and the horizontal x-axis can be obtained bywhere is the center angle of base thickness. The length of generating lines on points and can be obtained according to the gear geometry. Then the node on the theoretical tooth profile can be gotten by modifying the coordinate value of the node as the theoretical length of the generating line.
(a)
(b)
Since the position of the theoretical point is obtained by adjusting the length of the generated line, the mapping method of nodes on the flank can be conveniently applied to the tooth contact analysis of gears with profile modification.
5. Helical Gear Contact Analysis and Discussion
5.1. Boundary Conditions and Contact Pairs on the Contact Surfaces
Based on the basic parameters shown in Table 1 and mesh grid density control parameters, namely, the number of divisions of the fillet and the involute curve, the number of elements in tooth width, face width, and foundation depth directions, the 3D helical gear FEM can be established automatically. The hexahedral meshes in the contact area are well refined according to the refinement level, and all the nodes on the surface are mapped to the theoretical position. The high-precision 3D FEM of the helical gear has been established.
All nodes on the inner hub radius of the driving gear have been constrained radially, free rotation allowed. The nodes on the inner hub are coupled with the master node to make sure all the constrained nodes have the same rotation deformation. All the nodes on the hub of the driven gear are completely constrained in all directions. The input torque is applied as tangential nodal forces on the nodes located at the inner hub of the driving gear.
Majority of the previous studies utilized the gap-element method to simulate the gear contact behavior and approximate the contact deformation. In the current model, advanced surface-surface nonlinear finite element contact method and augmented Lagrangian contact algorithm are applied in LTCA. All the nodes on the tooth surface are selected in the contact region, and the element type CONTA174 is applied to represent contact and sliding between contact surfaces and deformable surfaces.
5.2. Comparison of the Coarse-Mesh Model and Well-Refined Model
Comparison of the developed 3D helical gear FEM with and without contact regions refinement is presented in this section.
5.2.1. Root Bending Stress
The correlation between the number of elements on the fillet curve and the tooth root bending stress is shown in Figure 12. As is demonstrated in the graph, the root bending stress increases rapidly with the increase of the number of elements on the root fillet curve when the number of elements on the root fillet curve is less than 9 and flattens off at a level of around 55.3 MPa. Since the bending stress gradient is relatively smooth, the coarse-mesh model is fine enough to evaluate the root bending stress.
5.2.2. Contact Pattern Comparison
Figure 13 shows the contact stress distributions on the tooth surfaces of the coarse-mesh helical gear engaged at the pitch point. The FEM is built with the element size recommended by Coy, , and the number of elements along the involute curve iswhere is the effective working arc length of the involute tooth profile between points ② and ③ (Figure 3).
(a)
(b)
(c)
As can be seen from Figure 13, the contact ellipses are formed when the load is mainly borne by a single row of nodes and discretely distributed along the contact line. Since the element on the contact surface is not refined, this contact pattern also appeared in the literature [32]. This may cause stress concentration in the contact simulation. But when the load is shared by two nodes, the contact stress on the tooth surface is much smaller. For the input load T = 76.2 Nm, the contact stress on the path of the contact line is shown in Figure 14. The maximum and mean value of the contact stress are 546.61 MPa and 377.32 MPa, respectively. The peak-to-peak amplitude is about 350 MPa. According to ISO 6336-2 [33], the contact stress on the surface of the helical gear is estimated to be 548.5 MPa. Although the maximum contact stress is close to the theoretical calculations from the Hertz theory, there is drastic fluctuation along the line of contact, and the average value is about 31.2% lower than the result from ISO. It is hard to say which value should be chosen to indicate the contact status. Obviously, the maximum contact stress extracted in this case does not correctly reflect the contact stress of the real tooth surface. The result is most likely to be underestimated than the theoretical value if there are not enough nodes in the contact width direction despite the fact that the entire gear FEM is too large and time expensive to solve.
It is worth mentioning that due to the excessive mesh density of the entire model, the calculation becomes very time consuming and takes about 86 to 116 minutes, depending on the density of the grid. The higher the mesh density, the more difficult and time-consuming the calculation. Obviously, the computational time of entire refinement model is unacceptable.
By applying the contact region refinement technology proposed previously in this paper, the element size of the noncontact area is 0.52 mm and the element size near the contact area is refined to 0.05 mm, the contact pattern is greatly improved, and the results of contact simulation are more reasonable. The contact stress distribution on the surface of the refined model is given in Figures 15 and 16, which show an obvious effect and supremacy over the coarse-mesh model. The FEM of the helical gear with contact zone refinement has continuous contact stress distribution on the tooth surface. 3D contact stress distribution of the refined-mesh model is shown in Figure 16. The abscissa axis donates the contact width, which is about 280 μm in this case. The vertical plane (“X = 0”) where the coordinate origin is located is the theoretical surface of action, and the geometric contact line appears on this plane. It can be found that the contact stress distribution is asymmetric to the surface of action. The maximum contact stress line is not located exactly at the theoretical contact line due to the tooth bending deflection and contact slippage in the contact process. The contact stress on the geometric contact line and actual contact line, as well as the sliding distance, can be figured out through measuring the contact stress on the surface as shown in Figure 17. The maximum and average contact stress on the actual contact line are about 573.3 MPa and 541.6 MPa, which is only about 1.26% lower than that of the theoretical calculations of the mean value. What stands out in this figure is the difference between the contact stress of the coarse-mesh model (black dashed line) and the well-refined model (red and blue solid lines). The upward trend of the graph is related to the change of the radius curvature on the contact surface. The contact edge relaxation effect can obviously be observed from the tooth contact analysis of the well-refined-mesh model.
As a comparison, the computational time of the local refinement model is much less with respect to the entire refinement model, which is reduced to about 20 minutes.
5.2.3. Grid Density and Normal Contact Stiffness Factor
The angular rotation of the driving gear can be obtained by retrieving the tangential displacement value of the gear hub (UYG), and then the transmission error is represented by the angle
The torsional mesh stiffness is given by
And the normal mesh stiffness along the line of action is given by
The grid density and normal contact stiffness have significant effects on the results of finite element contact simulation. The contact stress and mesh stiffness of the coarse-mesh finite element model with different mesh density and normal contact stiffness factor FKN are shown in Figures 18 and 19. Here, we can see that with the increase of the mesh density and the normal contact stiffness factor, the contact stress and the mesh stiffness increase and the mesh stiffness flattens off at a maximum value of around , while the contact stress is still increasing when the number of elements in the height direction increases to 45. Figure 20 shows that the computing time increases drastically with the increase in the number of elements on the involute but has little relationship with the FKN factor.
The contact stress and mesh stiffness of the mesh-refined model in the contact zone with 15 elements in the tooth height direction are shown in Figures 21 and 22. The refinement level is from 2 to 6, and the normal contact stiffness factor is from 0.01 to 4. From the graph, we can see that the maximum contact stress increases with the refinement level and FKN, while the mesh stiffness increases with the refinement level and decreases with FKN. Both the contact stiffness and the mesh stiffness reach a plateau eventually if the element size in the contact region is small enough and the FKN factor is large enough. The maximum contact stress and the mesh stiffness are 549.5 MPa and , respectively, which are in accordance with the theoretical calculation results. However, the computing time spent on the simulation is about 18 minutes at most, which is approximately 1/5 of the time demand in coarse-mesh contact model simulation.
Obviously, the local hexahedral mesh refinement in the contact zone can greatly reduce the simulation time and improve the accuracy of the calculation results.
5.3. 3D LTCA for the Helical Gear Pair
5.3.1. Contact Stress and Bending Stress Evaluation
Loaded tooth contact analysis and tooth stress calculation are conducted under a torque of 76.2 Nm for the gear pair with no misalignment error. The first solution is carried out at the mesh position of the pitch point (0°, as shown in Figure 23), then the gear is rotated with an angle increment, and the contact region refinement function is called to identify and refine the element in the contact area for the next solution. Simulations are carried out by using a looping program written by using MATLAB to calculate at different positions in the entire meshing cycle.
The bending stress and tooth contact stress of the refined-mesh model (indicated by the solid line) as well as the coarse-mesh model (indicated by the dashed line) of the whole mesh cycle are shown in Figures 23 and 24. According to ISO 6336-2, the root bending stress and contact stress on the surface of the helical gear are estimated to be 51 MPa and 538.5 MPa, respectively. In refined-mesh model, the maximum value of the bending stress and contact stress are about 60 MPa and 650 MPa, eliminating the effects of edge contact in alternating areas of single and double teeth in the mesh cycle, and the average value of bending stress and contact stress are about 46.2 MPa and 546 MPa. There are the maximum −9.4% and 1.39% difference between the ISO and the FEM results of the mean value. It can be observed a good agreement between the FE model and the analytical value.
While in the coarse-mesh model, the average value of maximum bending stress and maximum contact stress are about 46.06 MPa and 451 MPa, which are consistent with the observation from 5.2.1. The bending stress calculation has lower requirements on mesh density, and on the contrary, contact stress calculation is sensitive to the element size. The results of the contact stress show great difference as long as the mesh grid on the surface is not fine enough. It can be concluded that mesh density is critical for contact analysis.
It is worth noting that, in the alternating areas of single and double teeth, both the root bending stress and the tooth surface contact stress show anomalous values due to the edge contact effect.
5.3.2. Transmission Error and Tooth Mesh Stiffness
The static transmission error (STE) and mesh stiffness of individual tooth along with three teeth engagements of the refined-mesh and coarse-mesh helical gear model are calculated. STE of the external gear are shown in Figure 25, and the corresponding mesh stiffness is shown in Figure 26. The results of the refined-mesh model and coarse-mesh model are indicated by solid and dashed lines, respectively. The green lines in these figures are combined STE and mesh stiffness, and the others are individual tooth STE and mesh stiffness. It can be found from Figure 26 that the combined mesh stiffness is much smaller than the sum of the individual tooth mesh stiffness. Due to the nonlinear nature of contact deformation, the contact area increases as the load increases, and the contact stiffness increases as well. Therefore, the traditional method has far too high value in the prediction of combined mesh stiffness by parallel superposition of single-tooth mesh stiffness. Deformation coordination conditions for multitooth meshing must be taken into consideration. As can be seen from Figures 25 and 26, the STE of the coarse-mesh model is slightly larger than that of the refined-mesh model and the mesh stiffness is in the opposite way. This is because the contact deformation of the coarse-mesh model is larger than that of the refined-mesh model and the contact deformation accounts for a small proportion to the overall deformation. Therefore, the mesh density has little effect on the transmission error and the meshing stiffness.
5.3.3. Tooth Load Sharing Characteristics
By extracting the bearing force for each individual gear tooth at different meshing positions, quasistatic load sharing characteristics for each gear is determined by dividing the total normal load of every contact tooth with the sum of the normal load contributed by each pair in contact. Figure 27 shows the tooth load-sharing rate of the external gear pairs in the multimesh period. The results of the refined model and the coarse-mesh model are indicated by the solid and dashed lines, respectively. As can be seen from Figure 27, the load distribution curves of the two models are basically the same, except for the slight difference in the time of contact. The coarse-mesh model enters contact ahead of time. This is also because that the contact deformation of the coarse mesh is slightly larger than the actual value. The curve of the refined model is smoother and more accurate than that of the coarse-mesh model.
6. Conclusions
The method described in the present paper makes easy establishment of a 3D FEM of external helical gear automatically, high contact ratio gear included. The modelling process is parameterized and automated. Only the position parameter, the gear geometry parameters, and the FEM mesh density control parameters are needed in the program. Assembly errors such as center distance error and direction error can also be introduced into the modelling coordinate systems.
A high-precision mathematical model of tooth surface is the foundation to improve the accuracy of LTCA. The mathematical representations of the tooth profile including the root transit curve and the involute tooth flank of the external gear at the end face are obtained based on the manufacturing process and the equation of the cutter.
A bottom-up approach is used to generate the FEM: discrete nodes are created firstly and then eight-node hexahedral elements are defined directly, without building the solid model or importing solid elements, i.e., nodes and geometric model, from CAD software into the finite element solver. All of the procedures have been done in finite elements software to simplify the software environment and avoid data loss during model delivery.
Another important contribution of the article is the hexahedral mesh refinement method on the contact surface of the 3D helical gear model. The node mapping technique together with the tooth profile generation method guarantees all the nodes on the tooth flank are in exact geometrical position. This method can provide a good prediction of the load distribution, contact pattern, and mesh stiffness while reducing calculation time. The method presented in this paper makes it possible to get high-accuracy LTCA results on a personal computer. Therefore, it would have more hopeful prospects in gear loaded contact analysis.
Data Availability
All data used to support the findings of this study are available from the corresponding author upon request.
Conflicts of Interest
The authors declare that they have no conflicts of interest.
Acknowledgments
The authors gratefully acknowledge the support of the National Natural Science Foundation of China (no. 51505100), Natural Science Foundation of Heilongjiang for Young Scholars (no. QC2015046), and National Supercomputer Center in Shenzhen.