# Parameterized High-Precision Finite Element Modelling Method of 3D Helical Gears with Contact Zone Refinement

Authors: Yanping Liu; Yongqiang Zhao; Ming Liu; Xiaoyu Sun
Journal: Shock and Vibration
 (2019)

---

## 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.


---

## Body

## 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

Figure1(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 rct (part II) connects part I and part III. In the rigidly connected coordinate system σcoc−xc,yc of the rack cutter, the fillet curve is represented by the following equations:(1)xc=a−rctsinθ,yc=−b−rctcosθ,0<θ<π2−αt,where(2)a=πmt4+hat∗mttanαt+rccosαt,b=hat∗mt−rcsinαt−xtmt,rct=ct∗mt1−sinαt,where mt is the transverse module, αt is the transverse pressure angle, hat∗ and ct∗ are the transverse addendum coefficient and tip clearance coefficient, and xt is the transverse modification coefficient.Figure 1
Illustration of external gear generation. (a) Geometry parameters of a rack cutter in transverse sections. (b) Centrodes of the rack cutter and spur gear.
(a)
(b)Part III is represented by the following equation:(3)yc=−1tanαtxc+1tanαt⋅πmt4+xtmt.The superscriptt indicates that the values are measured in the transverse section. These equations are equally applicable to the spur gear.The coordinate systemsσcoc−xc,yc, σgog−xg,yg, and σ0o0−x0,y0 shown in Figure 1(b) are used to derive the tooth profile equations. The coordinate systems σc and σg 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 σc to σg is based on the following matrix equation:(4)xgyg1=Mgcxcyc1,and the translational matrix Mgc of the rack gear is(5)Mgc=cosφgsinφgrgsinφg−φgcosφg−sinφgcosφgrgcosφg+φgsinφg001.The generated profile of the gear in the coordinateσg 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 l of the rack cutter and the point xc,yc on the rack cutter have the relationship as(6)l=xc+yctanγ.In the fillet curve part (arc II) of the rack cutter, the tangent angleγ of the arc and its parameter θ are related by(7)tanγ=dycθdxcθ=−tanθ,where γ is the angle between the tangent line of the arc at point xc,yc with the x-axis.Substituting (1), (2), and (7) in (6) yields the following equation of the contact line represented by the arc parameter θ:(8)l=a+btanθ.According to Figure1, the rotation angle φg of the gear and the displacement l of the rack cutter are related by l=rgφg, and then the meshing equation of the rack cutter and the gear is represented by(9)rgφg−a−btanθ=0.Then, the equation of the transit curve of the gear represented by the arc parameterθ and gear rotation angle φg is(10)xg=acosφg−bsinφg−rcsinθ+φg−rgφgcosφg+rgsinφg,yg=−asinφg−bcosφg−rccosθ+φg+rgφgsinφg+rgcosφg,rgφg−a−btanθ=0.As to the oblique line (part III),(11)tanγ=dycdxc=−1tanαt.Substituting (3) and (11) into (6) yields(12)rgφg=xc−xtmttanαt+1tan2αtxc−πmt4.Then, the equation of the involute of gear represented by using the parameter of the rack cutter and rotation angleφg is given by(13)yc=−1tanαtxc+1tanαt⋅πmt4+xtmt,xg=xccosφg+ycsinφg+rgsinφg−φgcosφg,yg=−xcsinφg+yccosφg+rgcosφg+φgsinφg,rgφg=xc−xtmttanαt+1tan2αtxc−πmt4.


## 2.1. Transverse Tooth Profile of the External Gear

Figure1(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 rct (part II) connects part I and part III. In the rigidly connected coordinate system σcoc−xc,yc of the rack cutter, the fillet curve is represented by the following equations:(1)xc=a−rctsinθ,yc=−b−rctcosθ,0<θ<π2−αt,where(2)a=πmt4+hat∗mttanαt+rccosαt,b=hat∗mt−rcsinαt−xtmt,rct=ct∗mt1−sinαt,where mt is the transverse module, αt is the transverse pressure angle, hat∗ and ct∗ are the transverse addendum coefficient and tip clearance coefficient, and xt is the transverse modification coefficient.Figure 1
Illustration of external gear generation. (a) Geometry parameters of a rack cutter in transverse sections. (b) Centrodes of the rack cutter and spur gear.
(a)
(b)Part III is represented by the following equation:(3)yc=−1tanαtxc+1tanαt⋅πmt4+xtmt.The superscriptt indicates that the values are measured in the transverse section. These equations are equally applicable to the spur gear.The coordinate systemsσcoc−xc,yc, σgog−xg,yg, and σ0o0−x0,y0 shown in Figure 1(b) are used to derive the tooth profile equations. The coordinate systems σc and σg 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 σc to σg is based on the following matrix equation:(4)xgyg1=Mgcxcyc1,and the translational matrix Mgc of the rack gear is(5)Mgc=cosφgsinφgrgsinφg−φgcosφg−sinφgcosφgrgcosφg+φgsinφg001.The generated profile of the gear in the coordinateσg 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 l of the rack cutter and the point xc,yc on the rack cutter have the relationship as(6)l=xc+yctanγ.In the fillet curve part (arc II) of the rack cutter, the tangent angleγ of the arc and its parameter θ are related by(7)tanγ=dycθdxcθ=−tanθ,where γ is the angle between the tangent line of the arc at point xc,yc with the x-axis.Substituting (1), (2), and (7) in (6) yields the following equation of the contact line represented by the arc parameter θ:(8)l=a+btanθ.According to Figure1, the rotation angle φg of the gear and the displacement l of the rack cutter are related by l=rgφg, and then the meshing equation of the rack cutter and the gear is represented by(9)rgφg−a−btanθ=0.Then, the equation of the transit curve of the gear represented by the arc parameterθ and gear rotation angle φg is(10)xg=acosφg−bsinφg−rcsinθ+φg−rgφgcosφg+rgsinφg,yg=−asinφg−bcosφg−rccosθ+φg+rgφgsinφg+rgcosφg,rgφg−a−btanθ=0.As to the oblique line (part III),(11)tanγ=dycdxc=−1tanαt.Substituting (3) and (11) into (6) yields(12)rgφg=xc−xtmttanαt+1tan2αtxc−πmt4.Then, the equation of the involute of gear represented by using the parameter of the rack cutter and rotation angleφg is given by(13)yc=−1tanαtxc+1tanαt⋅πmt4+xtmt,xg=xccosφg+ycsinφg+rgsinφg−φgcosφg,yg=−xcsinφg+yccosφg+rgcosφg+φgsinφg,rgφg=xc−xtmttanαt+1tan2αtxc−πmt4.


## 3. Finite Element Modelling Methodology of 3D Helical Gear

The parameters of the gears are given in Table1.Table 1
Numerical data for the FE modelling of the gear pair.
Parameters
Gear
Pinion
Tooth number
29
23
Normal module (mm)
4
4
Normal pressure angle (°)
20
20
Helix angle (°)
23
23
Addendum coefficient
1
1
Dedendum coefficient
1.25
1.25
Face width (mm)
15
15
Shifting coefficient
0
0
Hub bore radius (mm)
35
25
Young’s modulusE (GPa)
206
206
Poisson’s ratioν
0.26
0.26The flowchart shown in Figure2 gives the parameterized finite element modelling method for the 3D external helical gear.Figure 2
Flow chart for creating the exact 3D helical gear FE model.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 Section2 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 A and B 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.Figure 3
Tooth profile division in the transverse section.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 ben1, on the involute contact profile between KP② and KP③ is n2, on the addendum circle between KP③ and KP⑨ is na, and on the web sideline between KP① and KP⑫ is nf. Keypoint KP① is set at the center rotation of the gear.According to the distribution of the KPs in Figure3 and the number of elements to be created on the corresponding lines, the number of KP①–KP⑮ is presented as(14)p1=2,p2=n1+2,p3=n1+n2+2+1=n1+n2+3,p4=n1+n2+2+1+1=n1+n2+4,p5=n1+n2+2+2+n1=2n1+n2+4,p6=2n1+n2+2+1=2n1+2n2+5,p7=2n1+n2+2+2=2n1+2n2+6,p8=2n1+n2+2+2+n1+1=3n1+2n2+7,p9=3n1+n2+2+1=3n1+3n2+7,p10=3n1+n2+2+2=3n1+3n2+8,p11=3n1+n2+2+2+nf=3n1+3n2+8+nf,p12=3n1+n2+2+2+nf+1=3n1+3n2+nf+9,p13=3n1+n2+2+2+2nf+1=3n1+3n2+2nf+9,p14=3n1+n2+2+2+2nf+2=3n1+3n2+2nf+10,p15=3n1+n2+2+2+3nf+2=3n1+3n2+3nf+10.Figure4(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]:(15)n2=sinvole,ebh=−0.2ce+1.2,0.9≤ce≤3,where sinvol is the arc length of the involute curve.(16)sinvol=rb2tan2αa−tan2α0,where c and e are the length and width of the element, respectively, and bh is the minimum Hertzian contact width in one mesh cycle:(17)bh=4Fnπbρ1ρ2ρ1+ρ21−υ12E1+1−υ22E2.Figure 4
Illustration of the finite element modelling process of the external helical gear. (a) KPs on the profile of one tooth. (b) KPs in the transverse plane. (c) Multilayer nodes of one tooth. (d) Single-tooth FEM. (e) Complete FEM of the driving gear. (f) Whole model of the 3D helical gear pair FEM.
(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 fillna 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(18)Nface=2n1+n2+nf⋅na+1.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 (Figure4(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 as(19)nb=we,where w is the face width of the gear.An hexahedral element is defined for every 8 neighbouring nodes to generate the FEM of one tooth (Figure4(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 Figure4(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σ18, 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 B1B2. The nodes Ni 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 NXNi=rb. 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., Ni−1 and Ni+1. 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.Figure 5
Coordinates parallel to the LOA.Figure 6
Illustration of the contact area and element group types: (a) contact strip; (b) the face of action; (c) points near the contact line are moved to the contact line and its sides; (d) standard contact region and the element groups.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 Figure6(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 Figure7. Once we get the minimum node number and give it to pc11, the number of other nodes in the remaining corner of the EG can be represented as(20)pc12=pc11+3n1+3n2+3nf+9,pc14=pc11+Nface,pc13=pc11+Nface+3n1+3n2+3nf+9,pci1=pci−11+1,pci2=pci−12+1,pci3=pci−13+1,pci4=pci−14+1,i=2,3,4.Figure 7
Three element group types. (a) Type I. (b) Type II. (c) Type III.
(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 Figure8. 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.Figure 8
Reconstruction method of type-III transition element group.
(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 Table2. 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 N filled on the edge along the contact line and the number of subelements En with the refinement level Lv are(21)N=1,2Lv−2,N=6,N2+3N+1,Lv=2,Lv>2.Table 2
Templates of hexahedral element refinement.
Refinement level
1
2
3
4
…
Refinement method
…
…
3D templates of hexahedral refinementThis refinement process is repeated along the contact line until all the original elements in the contact region are well refined, as is shown in Figure9. 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.Figure 9
Contact region refinement on the external gear flank.


### 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.Figure10 illustrates the coordinate system for node mapping. A cylindrical coordinate system σ0 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 σi can be gotten by translating σ0 along the z-axis of distance zi and rotating by angle ψi. Since xi-axis of the new system is still the symmetry axis of the tooth profile, the relationship between the rotation angle and the helix angle βn is(22)ψi=zitanβnr+i−12πz,where ei is the z coordinate value of the node on the tooth surface in σ0, r is the pitch radius, and i is the number of teeth from the reference position.Figure 10
Node mapping coordinate system.Pointp 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 p and theoretical point p″ of the external gear is shown in Figure 11(b). The angle θp between line Oip and the horizontal xi-axis can be obtained through the trigonometric function by extracting the polar radius ρp and vertical coordinate value NYp of the interpolation point p. Arc mapping point p′ is generated by the intersection of the arc with radius ρp and the theoretical involute. Since point p′ is on the theoretical tooth profile, the theoretical pressure angle, which has the relationship αp′=θp′+γ′, and θp′ can also be obtained by applying gear geometry theory. The offset angle between the arc mapping point p′ and point p is δ, δ=θp′−θp. The central angle φi of the node with a polar radius of ρi is calculated based on the theoretical tooth thickness formula of the gear. The angle between the generating line c′p′ and the horizontal x-axis can be obtained by(23)γ′=αp′+θp′−ψ2,where ψ is the center angle of base thickness. The length of generating lines on points p′ and p″ can be obtained according to the gear geometry. Then the node p″ on the theoretical tooth profile can be gotten by modifying the coordinate value of the node p as the theoretical length of the generating line.Figure 11
Mapping the node to the theoretical tooth profile.
(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.


## 4.1. Normalization of the Contact Region Elements

The coordinate systemσ18, 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 B1B2. The nodes Ni 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 NXNi=rb. 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., Ni−1 and Ni+1. 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.Figure 5
Coordinates parallel to the LOA.Figure 6
Illustration of the contact area and element group types: (a) contact strip; (b) the face of action; (c) points near the contact line are moved to the contact line and its sides; (d) standard contact region and the element groups.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 Figure6(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 Figure7. Once we get the minimum node number and give it to pc11, the number of other nodes in the remaining corner of the EG can be represented as(20)pc12=pc11+3n1+3n2+3nf+9,pc14=pc11+Nface,pc13=pc11+Nface+3n1+3n2+3nf+9,pci1=pci−11+1,pci2=pci−12+1,pci3=pci−13+1,pci4=pci−14+1,i=2,3,4.Figure 7
Three element group types. (a) Type I. (b) Type II. (c) Type III.
(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 Figure8. 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.Figure 8
Reconstruction method of type-III transition element group.
(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 Table2. 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 N filled on the edge along the contact line and the number of subelements En with the refinement level Lv are(21)N=1,2Lv−2,N=6,N2+3N+1,Lv=2,Lv>2.Table 2
Templates of hexahedral element refinement.
Refinement level
1
2
3
4
…
Refinement method
…
…
3D templates of hexahedral refinementThis refinement process is repeated along the contact line until all the original elements in the contact region are well refined, as is shown in Figure9. 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.Figure 9
Contact region refinement on the external gear flank.


## 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.Figure10 illustrates the coordinate system for node mapping. A cylindrical coordinate system σ0 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 σi can be gotten by translating σ0 along the z-axis of distance zi and rotating by angle ψi. Since xi-axis of the new system is still the symmetry axis of the tooth profile, the relationship between the rotation angle and the helix angle βn is(22)ψi=zitanβnr+i−12πz,where ei is the z coordinate value of the node on the tooth surface in σ0, r is the pitch radius, and i is the number of teeth from the reference position.Figure 10
Node mapping coordinate system.Pointp 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 p and theoretical point p″ of the external gear is shown in Figure 11(b). The angle θp between line Oip and the horizontal xi-axis can be obtained through the trigonometric function by extracting the polar radius ρp and vertical coordinate value NYp of the interpolation point p. Arc mapping point p′ is generated by the intersection of the arc with radius ρp and the theoretical involute. Since point p′ is on the theoretical tooth profile, the theoretical pressure angle, which has the relationship αp′=θp′+γ′, and θp′ can also be obtained by applying gear geometry theory. The offset angle between the arc mapping point p′ and point p is δ, δ=θp′−θp. The central angle φi of the node with a polar radius of ρi is calculated based on the theoretical tooth thickness formula of the gear. The angle between the generating line c′p′ and the horizontal x-axis can be obtained by(23)γ′=αp′+θp′−ψ2,where ψ is the center angle of base thickness. The length of generating lines on points p′ and p″ can be obtained according to the gear geometry. Then the node p″ on the theoretical tooth profile can be gotten by modifying the coordinate value of the node p as the theoretical length of the generating line.Figure 11
Mapping the node to the theoretical tooth profile.
(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 Table1 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 Figure12. 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.Figure 12
Tooth fillet bending stress.


#### 5.2.2. Contact Pattern Comparison

Figure13 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, e=bh=0.172mm, and the number of elements along the involute curve is(24)n2=sinvolbh=45,where sinvol is the effective working arc length of the involute tooth profile between points ② and ③ (Figure 3).Figure 13
Contact pattern of the coarse-mesh model.
(a)
(b)
(c)As can be seen from Figure13, 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.Figure 14
Contact stress on the contact line of the coarse-mesh model.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 Figures15 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.Figure 15
Contact stress distribution on the surface of the helical gear with contact region refined.Figure 16
3D contact stress distribution of refined helical gear model.Figure 17
Contact stress paths of the refined helical gear 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(25)θh=UYGrhub.The torsional mesh stiffness is given by(26)KT=Tθhub.And the normal mesh stiffness along the line of action is given by(27)KN=KTrb2.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 Figures18 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 1.45×105N/mm, 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.Figure 18
Contact stress of the coarse-mesh model.Figure 19
Mesh stiffness of the coarse-mesh model.Figure 20
Computing time of the coarse-mesh model.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 Figures21 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 1.53×105N/mm, 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.Figure 21
Contact stress of the refined model with 15 elements in the tooth height direction.Figure 22
Mesh stiffness of the refined model with 15 elements in the tooth height direction.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 Figure23), 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.Figure 23
Tooth root bending stress.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 Figures23 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.Figure 24
Tooth surface contact stress.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 from5.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 Figure25, 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.Figure 25
Transmission error of the external gear.Figure 26
Mesh stiffness of the external gear.


#### 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. Figure27 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.Figure 27
Tooth load sharing ratio of the external gear.


## 5.1. Boundary Conditions and Contact Pairs on the Contact Surfaces

Based on the basic parameters shown in Table1 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 Figure12. 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.Figure 12
Tooth fillet bending stress.


### 5.2.2. Contact Pattern Comparison

Figure13 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, e=bh=0.172mm, and the number of elements along the involute curve is(24)n2=sinvolbh=45,where sinvol is the effective working arc length of the involute tooth profile between points ② and ③ (Figure 3).Figure 13
Contact pattern of the coarse-mesh model.
(a)
(b)
(c)As can be seen from Figure13, 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.Figure 14
Contact stress on the contact line of the coarse-mesh model.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 Figures15 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.Figure 15
Contact stress distribution on the surface of the helical gear with contact region refined.Figure 16
3D contact stress distribution of refined helical gear model.Figure 17
Contact stress paths of the refined helical gear 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(25)θh=UYGrhub.The torsional mesh stiffness is given by(26)KT=Tθhub.And the normal mesh stiffness along the line of action is given by(27)KN=KTrb2.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 Figures18 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 1.45×105N/mm, 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.Figure 18
Contact stress of the coarse-mesh model.Figure 19
Mesh stiffness of the coarse-mesh model.Figure 20
Computing time of the coarse-mesh model.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 Figures21 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 1.53×105N/mm, 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.Figure 21
Contact stress of the refined model with 15 elements in the tooth height direction.Figure 22
Mesh stiffness of the refined model with 15 elements in the tooth height direction.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.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 Figure12. 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.Figure 12
Tooth fillet bending stress.


## 5.2.2. Contact Pattern Comparison

Figure13 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, e=bh=0.172mm, and the number of elements along the involute curve is(24)n2=sinvolbh=45,where sinvol is the effective working arc length of the involute tooth profile between points ② and ③ (Figure 3).Figure 13
Contact pattern of the coarse-mesh model.
(a)
(b)
(c)As can be seen from Figure13, 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.Figure 14
Contact stress on the contact line of the coarse-mesh model.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 Figures15 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.Figure 15
Contact stress distribution on the surface of the helical gear with contact region refined.Figure 16
3D contact stress distribution of refined helical gear model.Figure 17
Contact stress paths of the refined helical gear 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(25)θh=UYGrhub.The torsional mesh stiffness is given by(26)KT=Tθhub.And the normal mesh stiffness along the line of action is given by(27)KN=KTrb2.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 Figures18 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 1.45×105N/mm, 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.Figure 18
Contact stress of the coarse-mesh model.Figure 19
Mesh stiffness of the coarse-mesh model.Figure 20
Computing time of the coarse-mesh model.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 Figures21 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 1.53×105N/mm, 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.Figure 21
Contact stress of the refined model with 15 elements in the tooth height direction.Figure 22
Mesh stiffness of the refined model with 15 elements in the tooth height direction.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 Figure23), 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.Figure 23
Tooth root bending stress.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 Figures23 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.Figure 24
Tooth surface contact stress.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 from5.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 Figure25, 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.Figure 25
Transmission error of the external gear.Figure 26
Mesh stiffness of the external gear.


### 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. Figure27 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.Figure 27
Tooth load sharing ratio of the external gear.


## 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 Figure23), 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.Figure 23
Tooth root bending stress.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 Figures23 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.Figure 24
Tooth surface contact stress.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 from5.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 Figure25, 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.Figure 25
Transmission error of the external gear.Figure 26
Mesh stiffness of the external gear.


## 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. Figure27 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.Figure 27
Tooth load sharing ratio of the external gear.


## 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.



---

*Source: 5809164-2019-07-07.xml*